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ABSTRACT 


The  modeling  of  power  systems  has  been  primarily  driven  by  the  commercial  power 
utility  industry.  These  models  usually  involve  the  assumption  that  system  bus  voltage  and 
frequency  are  constant  However,  in  applications  such  as  shipboard  power  systems  this 
infinite  bus  assumption  is  not  vabd.  This  thesis  investigates  die  modeling  of  a  synchronous 
generator  and  various  loads  in  a  modular  fashion  on  a  finite  bus.  The  simulation  presented 
allows  the  interconnection  of  multiple  state-space  models  via  a  bus  voltage  modeL 

The  major  difficulty  encountered  in  building  a  model  which  computes  bus  voltage  at 
each  time  step  is  that  bus  voltage  is  a  function  of  current  and  current  derivative  terms. 
Bus  voltage  is  also  an  input  to  the  state  equations  which  produce  the  current  and  current 
derivatives.  This  creates  an  algebraic  loop  which  is  a  form  of  implicit  differential  equation. 

A  routine  has  been  developed  by  Linda  Petzold  of  Lawrence  Livermore  Laboratory 
for  solving  these  types  of  equations.  The  routine,  called  DASSL  (Differential/ Algebraic 
System  Solver),  has  been  implemented  in  a  pre-release  version  of  the  software  ACSL 
(Advanced  Continuous  Simulation  Language)  and  has  been  made  available  to  the  Naval 
Postgraduate  School  on  a  trial  basis.  An  isolated  power  system  is  modeled  using  this 
software  and  the  DASSL  routine.  The  system  response  to  several  dynamic  situations  is 
studied  and  the  results  are  presented. 
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L  INTRODUCTION 


A.  BACKGROUND 

Die  modeling  of  synchronous  generators  has  been  primarily  driven  by  the  commercial 
power  utility  industry.  These  simulations  most  frequently  make  the  assumption  that  the 
machine  to  be  modeled  is  connected  to  a  system  bus  in  which  voltage  magnimdr  and 
frequency  are  fixed  values.  This  so  called  "infinite  bus"  assumption  provides  good  results 
for  many  studies,  especially  those  involving  huge  power  grids.  However,  in  quite  a  few 
applications,  such  as  in  shipboard,  aircraft  and  isolated  emergency  power  systems  this 
assumption  is  not  valid. 

For  such  systems,  some  loads  may  be  a  significant  percentage  of  the  generator 
capacity.  When  a  large  load  is  started  on  an  isolated  generator,  neither  voltage  magnitude 
nor  frequency  remain  constant  In  these  situations  the  voltage  will  dip  appreciably  and 
possibly  cause  other  sensitive  loads  on  the  system  bus  to  fail. 

It  is  therefore  important  to  be  able  to  model  accurately  the  behavior  of  a  an  isolated 
power  system.  The  design  engineer  interested  in  building  the  smallest  least  expensive 
machine  that  will  do  the  job  needs  to  know  how  the  entire  system  will  behave  during 
dynamic  loading. 

B.  THESIS  OVERVIEW 

The  most  common  methods  for  studying  power  systems  and  the  interaction  between 
sources  and  loads  involves  use  of  the  infinite  bus  model.  The  interaction  between  sources 
and  loads  is  done  by  load  (power)  flow  analysis.  When  a  load  is  changed  in  such  a 
simulation  the  power  demand  must  be  satisfied  by  increasing  power  supplied  by  a  source. 
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This  approach  is  fine,  however,  it  still  assumes  each  independent  submodel  is  on  an  infinite 
bus.  So  for  example,  fluctuations  in  voltage  magnitude  and  frequency  occurring  in  (me 
generator  are  considered  negligible  in  the  overall  system. 

The  primary  goal  of  this  study  is  to  develop  a  system  model  which  allows  sources  and 
loads  to  be  connected  to  a  bus  voltage  model  which  accurately  reflects  the  bus  voltage 
behavior  during  transients  and  the  effects  of  the  transient  bus  voltage  on  the  loads  and 
sources.  The  approach  taken  is  to  develop  an  overall  system  model  which  consists  of 
accepted  accurate  stand-alone  source  and  load  submodels.  These  submodels  are  then  tied 
together  by  a  bus  voltage  model.  In  final  form,  the  system  model  allows  the  simple 
connection  of  multiple  sources  and  loads. 

Figure  1  represents  the  type  of  isolated  power  system  studied  in  this  thesis.  The 
system  consists  of  a  gas  turbine  prime  mover,  synchronous  generator,  a  system  bus  and 
system  loads.  Models  for  each  element  of  the  isolated  system  are  presented  independently 
then  all  are  tied  together  with  closed  loop  control  to  produce  the  system  model. 

1.  Synchronous  Generator  Model 

Much  work  has  been  done  to  develop  accurate  models  for  rotating  electrical 
machinery.  Krause  [Ref.  l,pp  211-267]  develops  a  state-space  model  for  a  synchronous 
machine  using  the  well  known  Park's  transformation  [Ref.  2].  Circuit  voltage  equations 
are  first  developed  in  the  three-phase  a-b-c  reference  frame.  The  subsequent  application  of 
Park's  transformation  has  the  advantage  of  referencing  all  state  variables  to  an  orthogonal 
(q-d-O)  reference  frame.  The  reference  frame  transformation  changes  time  varying  winding 
inductance  values  into  constants.  The  state-space  model  may  be  developed  with  either 
winding  current  or  magnetic  flux  linkage  as  states.  The  model  used  in  this  presentation  is 
formulated  with  current  states. 
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Figure  1.  Isolated  power  system. 
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2.  Load  Models 

Two  load  models  are  presented,  a  three-phase  resistive-inductive  (R-L)  load  and 
an  induction  motor  load.  These  models  are  developed  in  a  similar  manner  to  the 
synchronous  generator  model.  The  equations  describing  load  circuit  behavior  are 
transformed  into  the  same  orthogonal  reference  frame  as  that  used  for  the  synchronous 
generator.  Tire  load  model  state  equations  are  presented  with  current  states. 

3.  Bus  Voltage  Model 

By  formulating  the  submodel  state  equations  in  terms  of  current,  a  bus  voltage 
model  may  be  developed  based  on  satisfying  KirchofFs  Current  Law  (KCL)  at  the 
common  node.  Several  possible  models  for  the  bus  voltage  are  explored.  Ultimately,  a 
routine  for  solving  implicit  state  equations  and  algebraic  loops  will  be  introduced.  This 
routine  allows  a  bus  voltage  model  to  be  developed  which  supports  the  goals  of 
modularity  and  expandability. 

4.  Generator  Closed  Loop  Control 

Output  voltage  magnitude  and  frequency  must  be  controlled  so  that  bus  voltage 
will  remain  within  specification.  Voltage  control  is  accomplished  by  a  field  exciter  which 
senses  generator  output  terminal  voltage  and  adjusts  the  field  winding  excitation  voltage 
to  keep  terminal  voltage  at  the  desired  level  Frequency  control  is  accomplished  by 
driving  the  generator  with  a  prime  mover  which  is  under  the  control  of  a  speed  regulating 
governor.  Models  for  both  regulation  systems  are  presented. 

5.  Simulation  Software 

Speed,  ease  of  use,  quality  of  output  and  special  capabilities  were  considered 
when  choosing  the  simulation  software.  Work  was  done  in  the  programs  MATLAB  and 
SEMULINK  from  MathWorks  [Ref.  3]  and  in  ACSL  (Advanced  Continuous  Simulation 
Language)  from  Mitchell  and  Gauthier  Associates  [Ref.  4].  Both  are  excellent  for 
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modeling  systems  of  linear  and  non-linear  differential  equations.  However,  ACSL  was 
chosen  for  the  power  system  simulation  work  presented  here  due  to  the  special  capabilities 
of  this  package. 

A  pre-release  version  of  ACSL,  level  10F,  was  provided  to  the  Naval 
Postgraduate  School  on  a  trial  basis.  This  version  of  ACSL  contain*?  an  algorithm  for 
solving  differential  algebraic  equations  (DAEs)  which  is  described  in  Chapter  IV.  This 
algorithm  allows  systems  of  implicit  differential  equations  to  be  solved.  Specifically  of  use 
is  the  ability  to  solve  implicit  systems  formed  by  a  system  of  state  equations  subject  to  an 
algebraic  constraint  equation. 

5.  System  Model  Connection 

The  isolated  power  system  simulation  is  modular  in  concept.  It  consists  of 
several  submodels.  Each  submodel  is  a  stand-alone  model  which  is  tied  into  the  system  by 
the  bus  voltage  submodel.  The  source  and  load  models  are  well  understood  and  have  been 
validated  extensively  by  others.  The  bus  voltage  model  is  presented  and  validated  as  an 
independent  model  in  Chapter  IV.  After  each  piece  of  the  system  is  presented,  the  total 
isolated  power  system  is  developed  from  the  available  building  blocks.  The  ACSL  code 
for  the  system  model  is  described  in  detail 

6.  System  Model  Response  and  Validation 

After  the  system  is  connected  and  put  under  closed  loop  control,  it  is  validated 
by  comparing  it  with  a  finite  bus  system  model  developed  at  Purdue  University  by  Mayer 
and  Wasynczuk  [Ref.  5].  This  simulation  scenario  involves  starting  three  induction 
motors  on  a  system  bus  supplied  by  a  single  generator.  Plots  of  model  response  are 
presented  and  discussed. 

Additionally,  the  R-L  model  is  modified  for  the  case  where  the  resistive  part  of 
the  load  is  unbalanced.  The  system  model  is  exercised  by  operating  it  with  an  unbalanced 
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loading  condition.  Plots  of  the  model  response  to  this  condition  are  presented  and 
discussed. 

7.  Conclusions  and  Future  Work 

Finally,  conclusions  about  the  usefulness  of  die  finite  bus  model  are  presented. 
Suggestions  for  expanding  the  system  model  to  include  winding  saturation  effects,  more 
loads  and  a  parallel  generator  are  made.  The  need  for  more  effort  in  validating  the 
approach  is  also  discussed  along  with  some  suggestions  on  how  this  could  be 
accomplished. 
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H.  DEVELOPMENT  OF  THE  SYNCHRONOUS  GENERATOR  MODEL 


The  process  used  in  developing  the  model  for  a  synchronous  generator  is  as  follows. 
First  the  differential  equations  describing  the  circuit  behavior  of  each  winding  in  the 
machine  are  obtained.  Unfortunately,  because  both  the  current  and  inductance  terms  in 
the  equations  vary  with  time,  these  equations  are  complicated.  Using  Park's 
transformation  [Ref.  2],  the  equations  describing  the  machine  are  changed  to  an 
orthogonal  reference  frame  which  has  the  advantage  of  making  all  inductance  terms 
constant  This  transformation  from  machine  variables  to  reference-frame  variables,  along 
with  the  assumption  of  a  linear  relationship  between  current  and  flux  linkages,  allows  the 
model  state  equations  to  be  expressed  with  either  current  or  flux  linkages  as  the  states. 
For  a  synchronous  generator  the  orthogonal  reference  frame  used  will  be  the  rotor 
reference  frame. 

A.  SYNCHRONOUS  GENERATOR  EQUATIONS  IN  MACHINE  VARIABLES 
The  equations  used  to  develop  the  synchronous  generator  model  are  derived  from  the 
voltage  equations  for  the  windings  of  a  three-phase  machine.  These  equations,  which 
relate  voltage  to  current  and  magnetic  flux,  are  not  enough  to  completely  describe  the 
behavior  of  the  machine.  Additionally,  an  equation  is  needed  which  relates  rotor 
rotational  speed  to  torque  where  electrical  torque  is  described  as  a  function  of  current  or 
flux  linkage. 

1.  Machine  Variable  Voltage  Equations 

Figure  2  represents  a  two-pole,  three-phase,  salient-pole  synchronous  generator. 
The  as ,  bs,  and  cs  windings  are  on  the  stator  and  spaced  120°  apart  These  stator  windings 
are  identical,  sinusoidally  distributed  and  have  Nt  equivalent  turns.  On  the  rotor,  kq  and 
kd  are  damper  windings  while  th e/d  winding  is  used  for  applying  the  field  excitation.  The 
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rotor  windings  are  situated  in  an  orthogonal  q-d  reference  frame.  Rotor  windings  are  also 
sinusoidally  distributed  but  each  may  have  a  different  number  of  equivalent  turns;  N^,  Nu 
and  Np  respectively.  Note  that  the  current  direction  is  out  of  the  stator  windings 
(generator  convention).  The  series  resistive-inductive  pair  in  each  stator  and  rotor  circuit 
represents  the  electrical  characteristic  of  each  winding.  The  rotor  angular  position  and 
speed  are  represented  by  0rand  <or  respectively. 


bsaxis 


Figure  2.  Two-pole,  three-phase,  salient-pole  synchronous  machine. 
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By  summing  voltages  around  each  circuit  loop,  equations  (1)  through  (6)  are 
obtained.  Equations  (1),  (2)  and  (3)  represent  the  stator  windings  while  (4),  (3)  and  (6) 
model  the  rotor  circuits.  Voltage  in  the  inductive  elements  is  expressed  by  Faraday's  law 
where  the  induced  voltage  equals  the  rate  of  change  of  the  flux  linkages.  The  damper 
windings  are  short  circuited  at  the  ends  so  that  and  vu  equal  zero. 


v„  = 


v*.  = 


v_  = 


0  = 


= 


0  = 


f,,“  +  dt 


-r.it,  + 


0-L 

dt 
dt 

r,t°  dt 

dXu 

(fXt  fj 

rMiM  +  —r~ 
Mfi  dt 

Wu  + 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 


In  order  to  use  equations  (1)  through  (6)  to  develop  a  state-space  model  in 
terms  of  current,  the  flux  linkage  derivative  terms  must  be  looked  at  in  more  detail.  For  a 
linear  magnetic  system  the  flux  linkage  may  be  related  to  current  via  the  relationship 

X  =  Li  (7) 

where 

X  =  [X„  Xa  X^  X#  Xw]  (8) 

i  ~  [*«,  **»  ia  it,  ij*  <«]r  (9) 
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The  terms  of  die  matrix  L  represent  die  mutual  and  self  inductance  terms 
relating  flux  linkage  to  currents,  hi  general  would  relate  x-winding  flux  to  z-winding 
current.  So,  for  example,  the  expanded  expression  for  die  as-winding  flux  linkage  is 
written  as 


—  Lmi„  +  iw  +  Lmia  +  (10) 

where,  in  general,  both  inductance  and  current  are  functions  of  time.  Using  equation  (7 
the  inductance-current  product  may  be  substituted  for  flux  linkage  in  equations  (lj 
through  (6). 

Because  access  to  the  rotor  windings  is  difficult,  the  machine  parameters 
(winding  resistance,  inductances,  voltages  etc.)  are  most  commonly  referred  to  the  stator. 
Referring  values  from  the  rotor  to  the  stator  is  done  in  a  manner  similar  to  referring 
variables  from  the  primary  to  the  secondary  of  a  transformer  (via  the  turns  ratio).  Krause 
[Ref.  l:pp.  167]  uses  a  prime  to  denote  referred  variables.  In  this  derivation  all  variables 
may  be  assumed  to  be  referred  to  the  stator.  In  particular  the  inductance  matrix  which 
follows  as  part  of  the  compact  voltage  equations  is  written  in  referred  quantities. 

With  flux  linkage  expressed  in  terms  of  currents  and  winding  inductance,  a 
compact  vector-matrix  form  of  the  voltage  equations  may  now  be  written  as 


where  the  operator  p  represents  the  derivative  with  respect  to  time,  With  the 
equations  expressed  in  the  machine  reference  frame  the  derivative  operator  must  be 
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applied  to  the  inductance-current  product,  since  both  may  vary  with  time.  The  constant, 
is  a  function  of  referring  variables. 

In  equation  (11)  the  inductance  matrix  is  made  up  of  four  smaller  matrices.  The 
L,  matrix  relates  stator  winding  flux  linkage  to  stator  current  The  L,  matrix  relates  rotor 
winding  flux  linkage  to  rotor  current  The  matrix  relates  rotor  windings  to  stator 
windings.  These  matrices  have  the  following  form: 


L. 


Z*  +  La  -  L,  cos  20, 

cos  2(0,-|) 
-\la-L,  cos2(0,+|) 


-±Z*-Z,cos2(0,-|) 

A+A"A  00.2(0,-^) 

-It.- I,  cos  2(0,  +  *) 


cos  2(0,  +  |) 

-  4 «»2(0,  +  *) 
l*  +  LA  -  L,  a*  2(0,  + 

(12) 


i^cose, 

4*  cos(0,  -  ^) 

An,  C«5(®r  +  — ) 


sin  0, 

Aw  sin(0,  - 

Aw  sin(0,  +  y> 


sin  0, 

A-  sin(0f  - 

4w  sin(0,  +  ^) 


(13) 


Lr 


4*  +  4*  0  0 

0  Lyu  "*■  4w 
0  4w  Aw  +  Aw 


(14) 


The  inductance  terms  in  matrices  (12),  (13)  and  (14)  are  either  self  or  mutual 
inductances.  Elements  on  the  main  diagonal  of  the  L,  and  Lr  matrices  are  self  inductances 
and  are  made  up  of  leakage  and  magnetizing  inductance  parts  (Lk  +  A»).  Off  diagonal 
elements  of  L,  and  Lr  and  all  elements  of  L„  represent  mutual  inductances  and  therefore 
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are  assumed  to  have  no  leakage  inductance  part  The  variables  and  represent  the 
total  magnetizing  inductance  in  the  q  and  d  axes  respectively. 

Because  the  rotor  windings  are  wound  on  an  orthogonal  set  of  axes,  the  terms 
of  L,  are  easily  determined.  Orthogonal  magnetic  lines  of  flux  do  not  combine  so 
orthogonal  windings  have  zero  mutual  inductance.  The  mutual  inductance  for  windings 
which  share  a  common  axis  is  computed  as  for  a  transformer  which  is  die  product  of  die 
number  of  turns  divided  by  the  common  reluctance. 

The  situation  is  more  complicated  for  the  L,  and  L,,  matrices  The  inductance 
terms  depend  on  rotor  position.  Because  it  will  be  shown  dial  these  matrices  ate  greatly 
simplified  by  the  Park’s  transformation,  a  detailed  explanation  of  these  matrix  elements  will 
not  be  made  here  .  Figure  3  demonstrates  how,  in  the  case  of  a  salient-pole  machine, 
rotor  position  (which  changes  the  size  of  the  air  gap)  will  have  an  impact  on  the  reluctance 
path  and  therefore  on  the  inductance. 


■  magnetic  lines  efflux 


os- 


rotor  angle 


Figure  3.  Rotor  position  influence  on  winding  inductance. 
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For  the  as  winding,  minimum  reluctance  and  maximum  inductance  (LA  +  LB) 
occurs  when  9r  equals  90°  and  270°  while  maximum  reluctance  and  minimum  inductance 
(La  -  Lb)  is  experienced  at  0°  and  180°.  For  a  machine  with  a  round  rotor  LB  is  zero  and 
the  L,  terms  are  all  constant  A  complete  description  of  the  development  of  the  8r 
dependent  matrix  terms  may  be  found  in  Krause  [Ref.  l:pp  21 1-227]. 

2.  Machine  Variable  Torque  Equation 

The  torque  equation  is  developed  based  on  the  assumption  of  linearity  in  the  X  -i 
relationship  of  an  electromagnetic  system.  This  allows  the  energy  stored  (W)  in  the 
electromagnetic  system  to  be  expressed  by 

W  =  |  U2  (15) 


and  in  matrix-vector  form 


(16) 


Energy  or  work  is  also  the  product  of  force  and  displacement  Using  this  basic 
definition  yields 


W  =  TQ 
d W 

ae 


T  = 


(17) 


13 


where  T  is  torque  and  8  is  angular  displacement  Prom  equation  (17)  Krause  [Ref.  l:p 
217]  goes  on  to  My  develop  the  electrical  torque  equation  in  the  machine  reference  frame 
for  a  P- pole  generator 


d[L.  -  V] 
aef 


Laics  ^  (  Laics  ^ 


(18) 


The  electrical  torque,  is  positive  for  generator  action  when  the  stator  current  flows  out 
of  the  stator  terminals. 

One  more  differential  equation  is  needed  to  develop  a  state-space  model  Each 
voltage  equation  is  a  function  of  currents,  current  derivatives  and  rotor  position.  Rotor 
position  must  be  related  to  the  system  states.  The  second  derivative  of  rotor  position  may 
be  related  to  torque  which  in  turn  is  a  function  of  the  system  states.  The  final  differential 
equation  is  obtained  by  writing  the  relationship  for  the  mechanical  system  with  friction, 
windage  and  other  mechanical  losses  neglected 

iX»,  =  -  T.)  (19) 


where  7  is  the  inertia,  P  is  number  of  poles  and  T,  is  the  prime  mover  input  torque.  It  can 
be  seen  that  when  input  torque  is  greater  than  the  produced  electrical  torque  the  rotor 
acceleration  is  positive  and  the  machine  speeds  up.  A  large  load  will  cause  current  to  rise 
and  from  equation  (18)  electrical  torque  will  also  rise  resulting  in  deceleration  of  the 
machine. 
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B.  SYNCHRONOUS  GENERATOR  EQUATIONS  TRANSFORMED 


For  the  salient  pole  synchronous  generator  in  the  machine  reference  frame,  most 
inductance  terms  are  highly  dependent  on  0,  the  rotor  angle,  which  in  turn  is  dependent 
on  time.  This  makes  the  flux  linkage  derivative  term  a  complicated  chain  rule  expression 
(inductance  and  current  are  both  time  varying).  If,  however,  the  toms  of  the  inductance 
matrix  could  be  transformed  into  constants,  the  flux  linkage  derivative  could  be  simply 
expressed  as 


p\  =  L pi  (20) 

The  voltage  equations  may  be  rewritten  in  compact  form  by  assuming  that  such  a 
transformation  is  possible  and  making  the  substitution  suggested  by  equation  (20).  This 
yields 


y  as  ri  +  L pi  (21) 

Equation  (21)  then  becomes  the  basis  for  a  set  of  state  equations  describing  the  behavior 
of  the  synchronous  generator.  The  relatively  simple  form  of  equation  (21)  is  not  possible 
when  the  voltage  equations  are  expressed  in  the  machine  reference  frame,  in  other  words, 
with  stator  voltages,  currents  and  inductances  expressed  in  the  a-b-c  reference  frame. 

Fortunately  R.  H.  Park  developed  a  transformation  making  equation  (21)  possible. 
The  Park's  transformation  eliminates  time  varying  inductance  terms  and  introduces  a 
reference  frame  speed  term,  CO,  which  may  be  chosen  to  be  rotor  speed.  Thus  the  Park's 
equations  put  the  voltage  equations  in  the  orthogonal  q-d-0  reference  frame  of  the  rotor. 
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1.  Transformed  Voltage  Equations 

The  transformation  changes  variables  from  the  a-b-c  frame  to  the  q-d-0  frame  of 
reference.  For  an  arbitrary  vector  variable  /,  representing  voltage  or  current,  the 

transformation  matrix  K'  operates  as  follows 

u.  =  *'-U. 

where  the  transformation  matrix  is 


k; 


cos  0,  cos(6r  -  ~)  cos(0,  + 
sin  0,  sin(0r  -  sin(0r  + 

J 

1  1  1 


(23) 


and  rotor  position  is  defined  as 


er  =  £«>,($)* +  8,(0) 


(24) 


The  transformation  is  now  applied  to  the  voltage  equations  (12)  with  the 
following  result: 


Kr,pLt(K'yl 

|P(Lr)r(K;r‘ 


(25) 


where  the  transformation  applied  to  r,  yields 
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K;r.(K'rrl  =  K'r,I(K')"‘  =  r^KK)'1  =  rJL  =  r. 


(26) 


Because  it  has  equal  values  on  the  main  diagonal,  rg  is  not  altered. 

Multiplying  out  the  toms  of  the  inductance  matrix  gives  surprising  results.  By 
carefully  following  the  rales  of  matrix  multiplication,  applying  the  derivative  operator  m 
the  proper  sequence  and  using  die  correct  trigonometric  identities  the  voltage  equations 
may  now  be  expressed  as 


=  ~(r,  +  pLf)i„  -  1,00,1*  +  +  iwOM*  + 

V*  =  VM*  -  (r.  +  pLd)i*  -  KPAh  +  PKJm  +  PL*d'u 


v0,  =  -(»i  +  PK)hs 


0  =  -/’Wp  +  (r*  +  P4,)4, 


0  -  -pLmJds  +  PL*JjU  +  (rU  + 


(27) 

(28) 

(29) 

(30) 

(31) 

(32) 


where  all  inductance  terms  are  constants  and  die  following  definitions  apply 


L'  *  L*  +  ^ 

(33) 

Li  =  Lu  + 

(34) 

+ 

n 

(35) 

Am  ~  L#  +  Lm 

(36) 

Lu  =  L/h  + 

(37) 
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The  voltage  equations  may  be  used  in  the  form  of  equations  (27)  through  (32), 
but  it  is  more  common  to  see  the  inductances  expressed  as  reactances.  This  is 
accomplished  by  using  the  relationship  m hL  =  X.  The  inductance  is  multiplied  by  a  base 

angular  frequency  (often  60  Hz).  Machine  parameters  are  usually  provided  in  duns. 
Making  this  change  and  putting  die  state  equations  in  matrix  form  gives 


The  right  hand  side  of  equation  (38)  is  divided  into  three  parts 

v  =  A  Li  +  A  N(Ori  +  B pi  (39) 
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a  linear  part  (AJ),  a  nonlinear  part  (ANa>rfl  and  a  current  derivative  part  This  type  of 
equation  is  known  as  an  implicit  differential  equation.  This  is  because  a  particular  state 
equation,  for  example  equation  (27) 


v„  =  ~(r,  +  pLt)im  -  Ld(Ori *  +  pL^iH  +  iH  +  ^CDriw 


cannot  be  written  explicitly.  That  is  with  one  state  derivative  isolated  on  the  left  and  a 
combination  of  states  only  on  the  right  This  is  a  drawback  to  this  development  since  most 
simulation  software  prefers  die  explicit  form  for  the  differential  equations  of  the  model. 

2.  Transformed  Torque  Equation 

The  complete  form  of  the  final  state  equation  takes  shape  in  the  transformed 
reference  frame.  The  dependence  on  angular  displacement  has  been  eliminated  and 
replaced  with  a  speed  term  0)r.  The  torque-acceleration  equation  will  provide  this  final 
equation.  Krause  [  Ref.  l:p  227]  substitutes  the  reference  frame  transformation  into 
equation  (18) 


r.  -  (f)(K;rW 


(40) 


and  after  considerable  work  arrives  at 


Now  when  equation  (19)  is  used  for  the  final  state  equation,  the  speed  derivative  can  be 
expressed  as  a  function  of  the  other  states. 

C.  CHOICE  OF  STATES 

Most  developments  of  the  Park's  equations  arrive  at  a  set  of  state  equations 
expressed  in  terms  of  magnetic  flux  linkage.  Krause  [Ref.  2:pp.  177],  Anderson  [Ref. 
6:pp.  85-88]  and  many  others  express  a  preference  for  the  flux  linkage  expressions  over 
the  current  expressions  because  drey  have  explicit  form.  The  implicit  set  of  equations 
requires  that  a  matrix  inversion  be  performed  or  that  some  sophisticated,  and  often  slow, 
routine  be  used  to  solve  the  problem.  The  matrix  inversion  often  involves  a  poorly 
conditioned  matrix  and  methods  such  as  LU  decomposition  may  introduce  significant 
error. 

Sources  and  loads  on  a  common  system  bus  do  not  share  flux  linkage  but  do  share 
currents  (and  therefore  current  derivatives).  Because  the  goal  here  is  to  develop  a  model 
which  allows  an  entire  power  system  to  be  built  up  in  a  modular  fashion,  the  system 
submodels  will  be  expressed  in  terms  of  current  Also,  intuitively,  bus  voltage  must  have 
some  functional  relationship  to  the  bus  current.  Solution  of  the  finite  bus  problem  relies 
on  choosing  to  express  system  state  equations  in  terms  of  current 

In  order  to  put  the  system  in  explicit  form,  equation  (39)  is  manipulated  to  isolate  the 
state  derivative  on  the  left  hand  side  of  die  equation 

Pi  =  V(-AL)i  +  V(-AN)<ori  +  Vy  (42) 

where  V  =  B*1.  This  involves  the  matrix  inversion  mentioned  above  and  therefore  the 
possibility  of  error. 
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In  order  to  minimize  the  possibility  of  error,  the  B  matrix  was  inverted  symbolically 
so  that  the  matrix  toms  of  equation  (42)  could  be  expressed  as  functions  of  the  given 
machine  parameters.  MATHCAD  4.0  [Ref.  7]  was  used  to  perform  die  matrix  invasion 
and  surprisingly  the  terms  were  not  terribly  unwieldy.  Appendix  A  contains  the 
MATHCAD  output  showing  the  symbolically  inverted  matrix.  These  results  were  then 
used  to  write  the  final  form  of  the  state  equations  for  the  simulation.  The  explicit  form  of 
the  system  state  equations  may  be  seen  in  ACSL  simulation  code  of  Appendix  B  and  are 
described  in  Chapter  VI. 
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m.  DEVELOPMENT  OF  LOAD  MODELS 


The  next  step  in  the  process  of  modeling  a  finite  bus  power  system  involves  modeling 
the  system  loads.  In  order  to  have  a  system  simulation  in  which  loads  and  sources  can  be 
put  together  in  a  modular  fashion,  the  load  model  equations  are  developed  with  current 
states. 

Two  load  models  are  looked  at,  a  simple  R-L  load  and  an  induction  motor.  The  final 
simulation  allows  for  either  type  load  to  be  connected  to  the  bus  alone  or  for  both  type 
loads  to  be  connected  in  parallel 

The  choice  of  load  model  was  motivated  by  the  fact  that  even  in  isolated  systems, 
such  as  a  shipboard  system,  the  system  load  can  be  looked  at  as  a  nearly  constant  power 
factor  load  most  of  the  lime.  An  R-L  load  model,  allowing  for  the  time  varying  of  its 
resistive  and  reactive  parts,  adequately  simulates  many  loading  conditions. 

The  induction  motor  is  a  very  common  large  load  onboard  ship.  Fire  pumps, 
hydraulic  pumps  and  large  ventilation  fans  are  some  of  the  uses  for  induction  motors.  For 
this  reason  an  induction  motor  model  was  also  chosen  for  a  system  load. 

A.  THE  R-L  LOAD  MODEL 

The  three-phase  R-L  load  is  represented  by  Figure  4.  The  diagram  may  represent  a 
balanced  or  unbalanced  system.  That  is  the  resistance  and  inductance  of  each  phase  may 
or  may  not  be  equal,  hi  practice  considerable  effort  is  made  to  balance  the  system  load, 
however,  power  systems  are  frequently  subjected  to  unbalanced  loading  conditions.  The 
model  will  be  developed  for  a  balanced  load.  Methods  to  handle  the  imbalanced  case  will 
be  discussed  in  Chapter  VL 
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Figure  4.  Three-phase  R-L  circuit 


The  mathematical  development  for  the  load  parallels  the  generator  model 
development  First  the  voltage  equations  are  written  down  as 

v„  =  +  pk* 

v*5  =  Vi*  +  PK 

+  PK 

or  in  matrix  form  as 

V*bcs  =  TlUcl  +  Lf PUcl  (46) 

where  the  balanced  resistance  matrix  is 


(43) 

(44) 

(45) 


r,  0  0 

0  rt  0 

0  0  r, 


(47) 


and  the  balanoed  indnctmoe  matrix  is 

'1^  0  O' 

L,  -  0  0  (48) 

[o  0  1^ 

If  the  mutual  inductance  between  phases  is  assumed  negligible,  the  inductance  matrix  has 
terms  only  on  the  mam  diagonal.  It  is  a  fairly  simple  matter  to  include  off-diagonal 
(mutual  inductance)  tem*  the  matrix  stOl  will  have  no  time  dependent  terms.  Note 
also  that  the  derivative  operator  in  (45)  applies  only  to  die  current  since  there  are  no  time 
varying  terms  in  L|. 

Next  the  load  must  be  convened  to  die  same  reference  frame  (q-d-O)  as  die  generator 
model.  The  same  transformation  used  in  equation  (25)  is  applied  to  the  R-L  load 
equations.  These  equations  may  now  be  written  as 

Vo.  =  1U  +  KltMW'U.]  <49> 

Since  both  the  transformation  matrix  and  current  vary  with  time  the  product  rule  for 
differentiation  yields 

Vo,  =•  KIii(K,T,Voi  +  (KJLtf>(K;)-,)Vo,  +  KW  *;)-'/>  Vo<  (30) 

Using  the  same  approach  as  (26),  the  first  and  third  term  on  the  right  hand  side  are  easily 
obtained.  The  second  term  illustrates  how  the  speed  term  was  introduced  in  equations 
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(27)  through  (32).  The  p  operator  is  applied  to  the  transformation  matrix  inverse  with  the 
following  result 


cos  0,  sin  0,  1 

-  sin  0,  cos  0,  0 

p(K)~'  =  p 

008(0,  -  8ill(0,  —  ““)  1 

- 

-«ln(0, -y)  008(0,  -y)  0 

cos(0, +y)  «*“(»,+ 1 

-  011(0,  +  y)  008(0,  +  y )  0 

(51) 


then 


k;l.cd. 


-  sin  0r  cos  0,  0 

-sin(0,  --y)  cos(0r  -  y )  0 

-  sin(0r  +  y )  cos(0r  +  y )  0 


=  L,mr 


0  10 

-10  0 
0  0  0 


(52) 


Finally,  after  substituting  reactance  for  inductance,  the  state  equations  for  die  load  in  q-d-0 
reference  may  be  written  in  expanded  form  as 


(53) 

(54) 

Vo,  =  riht  +  -j-  Phi 
to* 

(55) 

Equations  (53)  through  (55)  represent  an  explicit  form  of  the  load  model  state 
equations.  These  equations  will  easily  "plug  in"  to  the  final  system  model. 
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B.  THE  INDUCTION  MOTOR  LOAD  MODEL 

Unlike  the  ample  R-L  load,  the  induction  motor  model  does  not  result  in  a  set  of 
state  equations  which  are  explicit  in  current  states.  The  induction  motor  model 
development  is  somewhat  like  the  synchronous  machine  model  The  development  here 
will  not  be  as  detailed  as  the  synchronous  machine  model  was.  However,  key  differences 
in  the  equations  based  on  machine  geometry  and  operating  theory  will  be  explained  before 
the  final  set  of  equations  is  presented. 

1.  Voltage  Equation  Development 

Figure  5  represents  the  two-pole,  three-phase,  induction  motor  load.  The  stator 
windings  of  this  machine  are  identical  sinusoidally  distributed  and  spaced  120°  apart  As 
for  the  synchronous  machine,  these  windings  are  designated  as,  bs  and  cs  each  having  Ns 
equivalent  turns.  The  rotor  arrangement  however,  is  considerably  different  from  the 
synchronous  machine  configuration. 

For  the  model  development  rotor  windings  are  considered  to  be  identical 
sinusoidally  distributed  windings  spaced  120°  apart  on  the  rotor  with  Nr  equivalent  turns. 
The  rotor  windings  are  designated  or,  br  and  cr.  The  rotor  displacement  and  speed  are 
represented  by  and  to^  respectively.  The  rotor  windings  are  all  shorted,  although 

machines  are  available  which  allow  the  rotor  windings  to  be  excited  externally. 

The  assumptions  are  not  entirely  valid  because  many,  if  not  most,  induction 
motors  are  of  the  squirrel-cage  variety.  In  this  type  of  machine  the  rotor  windings  consist 
of  metal  bars  laid  into  the  rotor  which  are  shorted  at  the  ends.  Krause  points  out  that  in 
most  cases  uniformly  distributed  rotor  bars  are  adequately  described  by  the  sinusoidal 
assumption  [Ref.  l:p  167]. 

One  obvious  difference  in  the  induction  motor  geometry  from  that  of  the 
synchronous  machine  is  the  shape  of  the  rotor.  Because  the  rotor  is  round,  none  of  the 
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terms  of  the  stator- stator  inductance  matrix  depend  on  Qm.  The  only  rotor  position 
dependence  is  seen  in  the  inductance  terms  relating  stator  to  rotor  windings.  Note  that  the 
subscript  for  rotor  position  and  speed  is  nn  to  differentiate  it  from  the  synchronous 
machine  rotor  speed,  (0r. 


bs  axis 


Figure  5.  Two-pole,  three-phase,  induction  motor. 
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Performing  •  Kkcfcoffs  Voltage  Law  (KVL)  sum  around  each  loop  allows  the 
voltage  equations  to  be  written  in  compact  form  aa 


(56) 


where  all  variables  are  referred  to  the  stator  via  the  turns  ratio.  The  zero  elements  of  the 
voltage  vector  are  due  to  the  fact  that  the  rotor  windings  are  shorted  on  the  ends. 

The  inductance  matrix  is  made  up  of  smaller  matrices.  In  die  expressions  below, 
the  leakage  inductances  are  1^,  and  for  the  stator  and  rotor  windings  respectively.  With 
all  variables  referred  to  the  stator  windings,  the  only  magnetizing  inductance  term 
appearing  in  the  matrices  is  4»  (the  stator  magnetizing  inductance).  The  matrices  consist 
of 


r  iii 

4  +  4.  -^4.  -^4. 

-{4.  4  +  4.  --4. 

~4.  -^4.  4  +  4. 
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4  +  4.  ~4.  -^4. 

4.  4  +  4.  ~  4. 
--4.  ~\L"  4  +  4. 


(57) 


(58) 
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CO8  0,.  006(0,.  +  y)  008(0*  -  y) 

*  L*  COS(0*  ~  y)  008  0*  CO8(0*  +  y) 

CO8(0*  +  y )  008(0*  -  y )  COS  0* 


The  voltage  equations  must  be  transformed  to  the  orthogonal  reference  frame. 
Applying  the  transformation  matrix,  K',  results  in 

_  !>.  .TU1J  k;l.  (k:  Vo.  1 


and  after  performing  the  multiplication  and  making  the  reactance  for  inductance 
substitution,  the  equations,  in  the  separated  form  of  equation  (38)  become 


‘rM  0  0  0  0  Olfi,,' 
0  r,  0  0  0  0  i* 

0  0  r,  0  0  0 

0  0  0  rt  0  0  ir 

0  0  0  0  t,  0  i* 

0  0  0  0  0  rr  Iq, 


Ojrr 


0  U 
0  b. 
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where  the  following  definition*  «pply 


XM  »  (62) 

Xm  .«,<**+ I*)  (63) 

Xn  -tM^+f*)  (64) 

a>,  *  (m, -<»„.)  (65) 


It  is  interesting  to  note  that  the  speed  teem  developed  when  the  transformation  is 
applied  to  the  L„  matrix  (a>8),  is  the  difference  between  the  reference  frame  speed 
(a>r)and  the  motor's  rotor  qjeed(a>m).  The  reference  frame  speed  usually  chosen  for 
induction  motor  simulation  is  the  bus  voltage  electrical  angular  frequency.  For  the  system 
model  being  developed  here,  electrical  frequency  is  defined  by  the  rotor  speed  of  the 
synchronous  generator.  This  «peed  difference  term,  known  as  slip  speed,  is  basic  to  the 
operation  of  the  induction  motor.  The  rotor  currents  will  be  zero  at  steady  state  unless  the 
reference  speed  and  rotor  speed  are  unequal  Rotor  current  will  be  shown  necessary  to 
produce  torque  in  the  next  section. 

One  more  equation  is  needed  to  complete  the  state-space  model  As  was  the 
case  for  the  synchronous  machine  model  the  speed  term  is  related  to  die  other  states  via 
the  torque  equation. 

2.  Torque  Equation  Development 

Using  the  argument  based  on  energy  stored  in  an  electromagnetic  system  the 
electrical  torque  developed  by  an  induction  motor  can  be  shown  to  be 

t.  =  («> 
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for  a  P- pole  motor  [Ref.  l:pp  169*170].  Equation  (66)  does  not  involve  L,  or  Lr  because 
only  the  L„.  matrix  is  dependent  on  O^,.  This  equation  shows  that  torque  will  be  zero 

when  rotor  current  is  zero. 

Application  of  the  transformation  matrix  to  (66)  yields 

t.  =  <«7> 


which  in  terms  of  currents  may  be  expressed  as 

=  (68) 


with  Te  positive  for  motor  action.  With  an  expression  for  torque  as  a  function  of  current 
states  the  final  state  equation  may  now  be  written. 

The  friction  and  windage  losses  are  once  again  neglected  and  the  differential 
equation  for  the  mechanical  system  is  written  down.  In  the  case  of  a  motor,  electrical 
torque  will  accelerate  the  rotor  while  applied  load  torque  (T,)  will  slow  the  rotor  down. 
The  equation  describing  this  is 


P<»r 


(69) 


where  Jm  is  the  inertia  of  the  motor's  rotor  and  P  is  the  number  of  poles. 

3.  Explicit  Form  of  the  Induction  Motor  Model 

The  voltage  equations  (61)  may  be  manipulated  in  the  same  manner  as  those  for 
the  synchronous  machine  in  order  to  put  them  in  explicit  form.  MATHCAD  was  again 
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used  to  symbolically  invert  the  matrix  associated  with  the  state  derivatives.  The  result  of 
the  matrix  inversion  is  contained  in  Appendix  A  and  the  full  form  of  the  state  equation 
may  be  seen  in  the  ACSL  code  of  Appendix  B  and  are  described  in  Chapter  VL 
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IV.  THE  BUS  VOLTAGE  MODEL 


The  next  step  in  building  a  total  system  model  based  on  Figure  1  involves  connecting 
the  source  model  with  the  load  model  (or  load  models  in  parallel).  The  physical  entity 
which  joins  load  and  source  is  the  system  bus.  Figure  6  is  a  block  diagram  of  the  system 
to  be  modeled  without  the  field  excitation  or  speed  regulator  loops  closed. 
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Flgure  6.  Isolated  power  system  block  diagram. 

The  source  and  load  models  have  been  discussed  in  some  detail,  now  a  model  for  the 
system  bus  voltage  must  be  developed.  The  simplest  approach,  and  the  approach  most 
frequently  taken,  is  to  assume  the  bus  voltage  is  of  fixed  magnitude  and  frequency.  This  is 
the  so  called  infinite  bus  assumption. 

Another  method  of  modeling  the  bus  voltage  is  to  attempt  to  develop  a  dynamic 
mathematical  expression  for  bus  voltage.  The  difficulty  with  this  approach  is  that  the 
previously  mentioned  algebraic  loop  problem  must  be  avoided.  As  a  practical  matter,  the 
simulation  code  cannot  use  bus  voltage  to  solve  for  the  current  derivative  terms  if  bus 
voltage  is  described  as  a  function  of  those  current  derivatives. 


33 


Two  methods  are  presented  for  dealing  with  the  algebraic  loop.  The  first  uses 
algebraic  manipulation  to  eliminate  current  derivative  terms  from  die  equation  describing 
bus  voltage.  The  second  method  involves  treating  the  total  system  as  a  large  implicit 
model  by  using  the  DASSL  algorithm  [Ref.  8]. 

A.  INFINITE  BUS  MODEL 

Use  of  the  infinite  bus  model  greatly  simplifies  the  power  system  simulation. 
Anderson  [Ref.  2:p.  26]  notes,  "A  major  bus  of  a  power  system  of  a  very  large  capacity 
compared  to  the  rating  of  the  machine  under  consideration  is  approximately  an  infinite 
bus".  Simulations  of  this  type  have  been  done  extensively  and  are  well  understood.  The 
generator  and  load  models  presented  in  Chapter  n  and  m  have  been  validated  as  infinite 
biu.  models.  Park  [Ref.  2],  Krause  [Ref.  1],  Anderson  [Ref.  6]  and  many  others  have 
demonstrated  the  validity  of  these  models  with  both  flux  linkage  and  current  states. 
However,  for  reasons  previously  mentioned,  the  finite  bus  model  will  not  be  used  for  the 
isolated  power  system. 

B.  PARALLEL  LARGE  RESISTANCE  MODEL 

Another  method  of  modeling  bus  voltage  so  that  it  may  be  used  as  a  varying  input  to 
all  the  system  submodels  is  to  connect  a  large  parallel  resistance  on  the  bus.  Figure  7 
shows  this  conceptually.  The  use  of  a  very  large  resistance  allows  source  current  and  load 
current  to  be  approximately  equal.  The  bus  voltage  may  be  computed  using  Ohm's  law 
and  then  fed  back  as  an  input  to  the  load  and  source  models.  This  approach  eliminates 
voltage  dependence  on  the  current  derivative  but  does  not  accurately  model  the  system. 
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Additionally,  before  any  other  computational  errors  are  accounted  for,  some  accuracy 
is  sacrificed  because  a  small  current  is  bled  off  by  the  resistor.  While  small  current 
may  not  be  significant  in  many  applications,  this  solution  method  is  not  as  satisfying  as  the 
methods  which  follow. 
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Figure  7.  Parallel  resistor  bus  voltage  modeL 


C.  MATHEMATICAL  EXPRESSION  FOR  BUS  VOLTAGE 

One  satisfying  method  of  solving  the  bus  voltage  model  dilemma  is  to  develop  a 
mathematical  expression  based  on  well  understood  and  accepted  theory.  At  the  node 
connecting  load  and  source  the  sum  of  the  currents  is  zero  by  KirchofF s  Current  Law 
(KCL).  It  also  follows,  because  differentiation  is  a  linear  operation,  that  the  sum  of 
current  derivatives  is  zero. 

The  algebraic  loop  difficulty  is  introduced  when  current  derivative  terms  show  up  in 
the  mathematical  expression  for  bus  voltage.  This  is  because  bus  voltage  is  an  input  to  the 
equations  from  which  current  derivatives  are  computed.  Most  simulation  software  will  fail 
when  algebraic  loops  are  encountered.  This  problem  can  be  addressed  by  finding  a  way  to 
eliminate  the  current  derivatives  from  the  bus  voltage  equation. 
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In  order  to  demonstrate  this  approach,  a  simplified  system  will  be  examined.  Figure  8 
represents  a  circuit  containing  a  voltage  source,  V»  and  passive  elements  Lh  Lj,  R,  and  R2. 
A  single  current,  i„  flows  in  the  circuit  and  the  differential  equation  describing  the  system 
is 

V,  =  (I,  +  +  («,+«,)(, 

dt  (70) 

» 

Later  in  the  chapter  this  circuit  will  be  placed  under  closed  loop  control.  For  control 
system  work,  the  s-domain  transfer  function  form  of  the  system  equation  is 


Figure  8.  Simple  circuit  to  demonstrate  bus  voltage  equation  representation. 

Once  the  current  is  known  it  is  possible  to  compute  the  dynamic  voltage  behavior  at 
the  node.  This  voltage,  V„,  is  the  voltage  appearing  across  the  elements  Lq  and  R2.  The 
differential  equation  for  this  quantity  is 
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V,  =  L,&  +  R,i, 


(72) 


The  transfer  function  from  current  to  node  voltage  is 


V.(s) 

I,(s) 


L*s  +  R2 


(73) 


and  the  transfer  function  relating  source  voltage  to  node  voltage  may  now  be  written  as 
the  product  of  (71)  and  (73) 

Vn(s)  I,(s)  _  V„(s)  _  l^s  +  Ri 

I, (s)  Vf(s)  V,(s)  (I,  +  4)s  +  (*,  +  *>)  (74) 


The  transfer  function  representation  of  the  system  for  node  voltage  given  source 
voltage  has  no  delay  (the  order  of  the  numerator  is  not  smaller  than  the  order  of  the 
denominator).  *  This  is  the  way  the  algebraic  loop  phenomenon  manifests  itself  in  the  s- 
domain.  Node  voltage  can  not  be  properly  fed  back  to  contribute  to  the  source  voltage 
without  adding  a  controller  delay.  As  an  open  loop  transfer  function,  equation  (74)  may 
be  simulated  in  many  software  packages.  The  system  current  and  node  voltage  may  also 
be  computed  using  equations  (70)  and  (72). 

The  ultimate  goal  of  the  system  simulation  being  developed  here  is  to  be  able  to  take 
independent  submodels  and  tie  them  together  with  a  bus  voltage  model  This  modular 
approach  has  the  advantage  of  not  requiring  the  entire  model  to  be  redeveloped  if  one  of 
the  sources  or  loads  changes.  It  is  not  desirable  in  complicated  systems  to  develop  a  new 
system  model  each  time  a  submodel  changes,  figure  9  shows  how  the  baric  circuit  of 
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Figure  8  may  be  broken  m  two  pieces,  source  and  load,  to  demonstrate  die  modular 
approach. 


Li  Rj  Lg 


Figure  9.  Dividing  the  system  into  submodels. 


The  equations  describing  each  submodel  may  then  be  obtained  independently.  For 
the  source  submodel 


(75) 


These  submodels  are  accurate  and  it  is  easy  to  see  that  if  terminal  (bus)  voltage,  Vr  is 
a  fixed  value  (infinite  bus)  both  equations  are  easily  solved.  However,  a  solution  which 
models  the  source-load  interaction  and  the  terminal  voltage  dynamic  behavior  is  desired. 
In  order  to  do  this,  terminal  voltage  must  be  solved  at  each  time  step  of  the  simulation. 
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The  terminal  voltage  is  then  fed  without  delay  to  the  set  of  equations  describing  source 
and  load.  The  previously  derived  expression  for  node  voltage,  equation  (72),  can  not  be 
used.  A  current  derivative  term  appears  in  this  expression  which  will  cause  an  algebraic 
loop  error. 

One  possible  solution  is  to  develop  an  expression  which  eliminates  the  derivative 
terms.  Based  on  KCL  at  the  connecting  node  the  current  derivatives  may  be  set  equal 
The  resulting  equation  for  terminal  voltage  has  no  derivative  terms 


(77) 


This  allows  V,  to  be  solved  at  each  time  step  of  the  simulation  and  then  fed  into  the 
submodel  equations.  The  system  of  Figures  8  and  9  was  modeled  in  ACSL.  The  model 
was  formulated  using  both  the  single  system  approach  and  the  modular  system  approach. 
Appendix  B  contains  the  ACSL  code  used  for  this  simulation.  The  model  parameters  used 
were 


^  =  1.0ft  Rt  =  5.0ft 

4  =  0. 6H  Li  =  0. 2H 

and  the  source  voltage  was  a  .5  duty  cycle  square  wave  pulse  train  with  magnitude  equal 
to  1.0  and  a  period  of  2.5  seconds. 

Figure  10  compares  the  circuit  response  using  the  single  loop  model  with  the 
response  obtained  from  the  submodel  approach.  The  single  loop  solution  for  voltage, 
VNODE,  is  computed  from  equation  (72)  after  the  differential  equation  for  current  (70)  is 
solved.  The  plot  labeled  VT  is  the  solution  for  the  terminal  voltage  using  the 
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submodel  approach  and  equation  (77).  A  voltage  is  computed  at  each  time  step  then  fed 
as  an  input  to  fee  source  and  load  submodel  equations  (75)  rad  (76).  The  difference 
between  the  two  solution  methods  is  DELV. 

Comparison  of  fee  current  solution  in  each  submodel  is  also  significant  Ideally,  fee 
two  currents  will  be  equal  since  the  submodels  share  a  common  node.  DELI  is  the 
difference  between  i,  and  ij  and  DELID  is  the  variation  between  pix  and  pi?  Also  (dotted 
is  i,,  representing  current  flowing  in  the  system. 

Figure  10  demonstrates  that  the  bus  voltage  equation  submodel  yields  good  results. 
The  magnitude  of  error  in  voltage  and  current  is  extremely  small  and  this  approach 
achieves  the  desired  modularity  of  the  system  model.  There  are,  however,  some 
disadvantages  to  this  approach. 

The  current  error  variation  (DELI),  although  small,  is  still  growing  at  fee  end  of  fee 
run.  This  is  due  to  fee  formulation  of  fee  bus  voltage  equation.  The  current  derivatives 
are  set  equal  but  nothing  in  fee  equation  forces  fee  currents  to  stay  together. 

Another  problem  is  that  the  bus  voltage  equation  is  not  simple  for  complex  systems. 
The  equations  for  fee  synchronous  generator  wife  an  R-L  load  are  quite  complicated  and 
with  an  induction  motor  load  are  more  so.  This  approach  also  requires  that  the  bus 
voltage  model  be  redeveloped  for  different  source  and  load  submodels.  The  addition  of 
submodels  in  parallel  on  the  bus  creates  difficulties  for  fee  same  reasons  as  just  mentioned 
for  complex  systems. 

D.  DASSL  BUS  VOLTAGE  MODEL 

In  order  to  have  a  system  model  which  is  truly  modular,  a  better  solution  than  the  bus 
voltage  equation  must  be  found.  The  requirement  to  find  a  new  and  often  complicated 
expression  for  the  bus  voltage  each  time  a  submodel  is  altered  becomes  extremely 
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unwieldy  and  inconvenient.  The  previously  mentioned  DASSL  algorithm  provides  a  more 
satisfactory  solution  to  the  problem. 

1.  Row  DASSL  Works 

The  basic  idea  of  the  DASSL  routine  is  to  replace  the  derivative  in  die  implicit 
equation  with  a  difference  approximation  and  then  solve,  using  Newton's  method,  for  the 
derivative  at  die  current  time  step.  It  has  been  implemented  in  the  commercially  available 
simulation  package  ACSL.  The  class  of  problem  DASSL  is  designed  to  solve  are  known 
as  differential  algebraic  equations  (DAE).  DAEs  include  systems  of  implicit  differential 
equations  and  systems  of  equations  containing  algebraic  loops. 

To  show  how  die  routine  works  the  load-source- bus  problem  will  be  set  up  as  a 
DAE  or  implicit  equation.  One  way  to  write  the  general  DAE  [Ref.  9:p.  41]  which 
applies  to  the  case  of  die  load,  source  and  bus  voltage  models  is  to  represent  the  load  and 
source  model  equations  as  a  system  of  differential  equations  of  die  form 

Pi  =  /i(i.v.O  (78) 


coupled  to  an  algebraic  constraint 


Q  = 


(79) 


If  it  is  then  assumed  that  voltage  is  some  unknown  function  of  the  state  and  state 
derivative  vectors 


v  =  h(pi,  £,  r) 


(80) 
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the  system  equations  may  be  represented  in  a  hilly  implicit  fonn  by  substituting  for  voltage 
in  equations  (78)  and  (79).  Implicit  equations  of  this  type  have  been  solved  most 
successfully  using  backward  differentiation  formulas  (BDF)  [Ref.  9:p.42]. 

The  simplest  BDF  method  involves  replacing  the  derivative  with  a  first  order 
backwards  difference.  DASSL  extends  the  idea  of  die  BDF.  Rather  than  using  the  first 
order  difference,  the  derivative  is  approximated  by  the  *4  order  difference  and  k  ranges 
from  one  to  five.  Order  and  step  size  are  automatically  selected  based  on  the  behavior  of 
the  solution.  Because  of  the  flexibility  of  the  routine,  DASSL  has  been  shown  to  be  highly 
stable  and  robust  [Ref.  9:pp.  115-116]. 

Two  criteria  must  be  met  to  successfully  solve  a  DAE  system.  The  system  must 
be  solvable  and  implementable  [Ref.  9:p.  16].  Solvable  means  that  the  states  are 
differentiable  over  the  time  interval  of  interest  The  solution  must  be  smooth  enough  to 
make  this  possible.  Implementable  refers  to  the  solution  technique.  The  method  used  to 
solve  the  nonlinear  system  of  equations  must  provide  a  solution  at  each  time  step.  In  the 
case  at  hand,  as  is  often  done,  the  solvability  will  be  assumed.  The  current  states  are 
relatively  smooth  functions.  The  implementability  of  the  DASSL  routine  is  also  assumed. 

The  routine  as  implemented  in  ACSL  may  be  used  to  solve  implicit  differential 
equations,  algebraic  loops  and  systems  of  differential  equations  with  an  algebraic 
constraint  Conceptually  the  bus  voltage  problem  may  be  looked  at  as  a  system  of  state 
equations  with  an  algebraic  constraint  based  on  KiichofF s  Current  Law.  In  order  to  solve 
for  a  node  voltage  in  ACSL  using  the  implicit  equation  solver  the  constraint  equation  must 
first  be  defined.  Based  on  KCL  at  a  connecting  node,  a  residual,  r,  may  be  defined 

r  =  I  pi*  +  1 1*  (81) 

k  k 
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where  it  is  desired  to  keep  the  residual  at  or  near  zero.  The  node  voltage,  which  is 
assumed  implicitly  related  to  this  KCL  equation,  may  then  be  found  by  using  the  DASSL 
implicit  solver.  The  operator  in  ACSL  is  IMPLC,  and  the  syntax  is 

*  IMPLC(r,  (82) 

where  V-rrf r_fr  is  the  node  voltage  initial  condition  [Ref.  10:pp.  1-4]. 

2.  The  Advantage  of  the  DASSL  Bos  Voltage  Model 

Tire  fact  that  the  constraint  equation  (81)  for  the  system  is  simple  regardless  of 
the  complexity  of  subsystems  connected  at  the  node  of  interest  gives  the  DASSL  bus 
voltage  model  a  great  advantage  over  the  mathematical  model  presented  in  the  previous 
section.  It  is  a  simple  matter  to  add  or  remove  loads  and  sources  from  the  larger  system 
model  No  reformation  of  the  bus  voltage  model  is  required.  Additionally,  the  constraint 
equation  keeps  the  error  between  submodel  currents  close  to  zero.  Even  if  a  transient 
occurs  which  makes  the  error  grow  momentarily,  the  implicit  solver  adjusts  voltage  to 
move  the  current  error  back  towards  zero. 

The  simple  system  of  figure  9  was  simulated  using  the  implicit  feature  of  ACSL. 
The  same  model  parameters  were  used  as  in  the  simulation  results  presented  in  figure  10. 
The  code  used  may  be  seen  in  Appendix  B.  Figure  1 1  contains  the  simulation  results  using 
the  DASSL  routine.  As  with  the  results  of  the  bus  voltage  equation  submodel  the 
DASSL  bus  voltage  submodel  is  in  excellent  agreement  with  the  single  loop  solution. 
Notice  the  improvement  in  current  derivative  error  and  current  error,  DELID  and  DELI 
respectively. 
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Figure  11.  Comparison  ot  single  loop  solution  with  solution  obtained  using  the 
DASSL  bus  voltage  submodel  approach 


E.  THE  BUS  VOLTAGE  MODEL  AND  CLOSED  LOOP  CONTROL 

One  final  characteristic  required  of  the  bus  voltage  model  is  that  its  output  be  able  to 
be  fed  to  a  bus  voltage  control  circuit  The  power  system  model  which  is  being  developed 
here  requires  that  field  excitation  for  die  generator  be  controlled  in  order  to  maintain  the 
bus  voltage  at  or  near  specification.  To  accomplish  this,  the  bus  voltage  submodel  must 
accurately  track  the  terminal  voltage  behavior  during  transients  so  that  the  control  circuit 
submodel  responds  correctly. 

In  order  to  demonstrate  the  behavior  of  the  bus  voltage  submodel,  the  source-load 
system  of  Figure  9  will  again  be  used.  The  transfer  function  relating  source  voltage  to 
terminal  voltage  was  developed  and  given  as  equation  (74).  Using  this  transfer  function 
form  of  the  system  model  a  control  system  may  be  developed.  In  this  case  a  cascade 
controller  was  designed  using  the  root  locus  technique  which  produced  a  highly  damped 
response  with  a  fast  settling  time  and  small  steady  state  error.  The  transfer  function 
design  is  shown  in  Figure  12. 


Figure  12.  Transfer  function  form  of  the  simple  source- load  system  with  bus  voltage 

control. 


The  closed  loop  response  was  modeled  three  ways  using  ACSL.  The  transfer 
function  form  of  the  system  model  was  compared  with  a  closed  loop  form  using  both  the 
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bus  voltage  equation  and  the  DASSL  routine.  In  the  simulation  runs,  system  parameters 
were  initially  set  to  the  values  given  on  page  39.  The  voltage  reference  was  then  step 
changed  from  zero  to  one.  All  forms  of  the  system  model  produced  a  similar  step 
response.  Next  the  load  parameters  were  step  changed  to  R2  =  0.01  and  =  0.001.  The 
two  submodel  forms  of  the  system  exhibited  similar  behavior,  a  large  increase  in  load 
current  and  a  large  initial  dip  in  bus  voltage.  This  transient  was  followed  by  bus  voltage 
recovery  to  the  commanded  value. 

Figure  13  shows  the  response  of  the  closed  loop  system  to  a  step  change  in  the 
voltage  reference  from  zero  to  one.  Results  for  both  bus  voltage  submodels  are  shown. 
The  plot  labeled  DELV  is  the  difference  between  the  transfer  function  solution  and  the 
indicated  submodel.  Both  bus  voltage  submodels  provide  very  good  results  and  are  in 
excellent  agreement  with  the  transfer  function  model.  Note  the  plot  of  DELI  which  is 
moving  away  from  zero  in  the  bus  voltage  equation  submodel. 

Figure  14  is  the  response  of  the  system  to  the  step  load  change.  The  DASSL  results 
are  in  very  good  agreement  with  the  bus  voltage  equation  submodel  The  difference  to 
note  is  again  the  DELI  plot.  Figure  14  (b).  The  simulation  was  allowed  to  continue  to  100 
seconds.  The  DASSL  routine,  while  allowing  some  error  between  load  and  source 
current,  keeps  forcing  the  error  back  toward  zero.  The  bus  voltage  equation  allows  the 
DELI  error  to  grow  50  to  100  times  larger  than  the  DASSL  routine  allows. 

F.  CHOICE  OF  BUS  VOLTAGE  MODEL 

While  the  methods  presented  here  are  by  no  means  exhaustive,  several  reasonable 
possibilities  for  a  bus  voltage  model  have  been  looked  at  Some  advantages  and 
disadvantages  of  each  have  been  discussed.  Remember  that  modularity  and  simplicity  are 
desired  characteristics  of  the  total  system  model 
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Figure  13.  System  response  of  both  voltage  submodels  to  step  change  in  reference 
voltage.  DELV  compares  model  response  to  a  transfer  function  form  of  die  system. 
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Figure  14.  Step  change  in  load  for  both  submodels;  (a)  terminal  voltage  and  currant 
response,  (b)  source  and  load  current  difference  to  100  seconds. 


The  only  choice  presented  which  allows  both  modularity  and  simplicity  in 
implementation  as  the  total  system  model  grows  is  the  DASSL  bus  voltage  model  The 
routine  accurately  solves  for  bus  voltage  and  then  feeds  that  solution  to  all  connected 
submodels.  Additional  submodels  may  be  included  by  simply  adding  the  current  and 
current  derivative  terms  from  the  new  submodel  into  the  KCL  constraint  equation.  The 
DASSL  routine  is  easy  to  use,  fast,  accurate  and  robust 
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V.  SYSTEM  MODEL  DESCRIPTION 


Most  of  the  pieces  required  to  build  an  isolated  power  system  operating  on  a  finite 
bus  have  been  developed  and  presented.  Validated  submodels  for  source,  loads  and  bus 
voltage  are  available  to  be  used  as  building  blocks  for  a  larger  modeL  This  section  puts 
the  system  model  together  by  providing  a  means  for  closed  loop  control,  a  description  of 
the  per-unit  (pu)  system  and  a  detailed  description  of  the  ACSL  code  used  in  the  final 
system  simulation. 

A.  SYSTEM  CLOSED  LOOP  CONTROL 

In  order  to  complete  the  system  model  and  put  it  under  closed  loop  control,  two 
more  submodels  must  be  developed.  These  are  the  field  excitation  system  and  the  prime 
mover  and  speed  governor  system.  The  field  excitation  system  is  a  voltage  regulator  and 
exciter  which  senses  the  magnitude  of  the  bus  voltage  and  then  adjusts  the  voltage  applied 
to  the  field  winding  as  necessary  to  maintain  the  bus  voltage  at  or  near  the  commanded 
reference  level.  The  speed  governor  senses  the  mechanical  speed  of  the  rotor  and  applies 
a  control  signal  to  the  prime  mover  which  maintains  rotor  speed  at  or  near  the  commanded 
reference  level. 

1.  Field  Excitation  System  Model 

There  are  many  types  of  field  excitation  systems  used  for  synchronous  generator 
voltage  control.  The  IEEE  Type  2  representation  [Ref.  10]  presented  here  is  typical  for 
power  systems  used  aboard  US  Navy  ships.  It  is  also  the  type  of  field  excitation  system  in 
the  model  which  is  used  for  comparison  in  the  next  chapter. 

Figure  IS  is  a  block  diagram  of  the  IEEE  Type  2  regulator  and  exciter.  The  first 
section,  from  the  summing  junction  to  Vn,  is  the  regulator.  From  the  second  summing 
junction  to  the  output  is  the  exciter  section.  There  are  two  saturation  functions,  one 
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associated  with  each  part  of  the  system.  The  details  of  implementing  this  model  in  the 
simulation  are  discussed  in  the  section  describing  the  ACSL  code. 


Figure  15.  IEEE  Type  2  regulator  and  exdter  system  based  of  Ret  10. 

2.  Prime  Mover  and  Speed  Governor  Model 

As  for  the  excitation  system,  many  models  are  available  for  prime  mover  and 
speed  governor.  Steam  turbines,  hydro- turbines,  diesel  engines  and  gas  turbines  are  a  few 
examples  of  the  types  of  system  used  to  drive  synchronous  generators.  A  simple  transfer 
function  form  for  one  of  these  systems  might  consist  of  two  simple  first  order  delays,  one 
for  the  speed  governor  and  one  for  the  prime  mover.  Choice  of  the  prime  mover  and 
governor  model  was  driven  by  the  desire  to  compare  results  of  this  system  simulation  with 
other  work 

Mayer  and  Wasynczuk  [Ref.  5]  provide  a  model  for  an  Allison  501  gas  turbine 
engine  and  governor.  Figure  16  is  the  s-domain  representation  of  this  model  The  model 
is  a  per-unit  model  A  per-unit  model  references  all  model  variables  to  some  base  value. 
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Per-unit  models  are  very  common  in  power  system  iumahtion  work  and  wfll  be  described 
in  more  detail  in  the  next  section. 


Figure  16.  Gas  turbine  and  speed  governor  model  based  on  Ref.  5. 

B.  THE  PER-UNIT  SYSTEM 

The  per-unit  (pu)  system  was  developed  to  simplify  the  calculation  and  interpretation 
of  results  in  power  system  simulation  work.  A  pu  quantity  is  defined  as  an  actual  quantity 
(voltage,  current,  power  etc.)  divided  by  a  base  or  reference  quantity.  The  base  values  are 
selected  according  to  known  characteristics  of  the  machine  being  studied.  The  final 
system  simulation  presented  here  is  in  the  per-unit  system. 

The  pu  system  is  based  on  machine  rated  bus  voltage,  power  and  synchronous  speed. 
Machine  parameters  are  usually  supplied  by  the  manufacturer  in  pu.  The  base  quantities 
used  here  are 

W  rated  operating  voltage 
PbaM  :  rated  volt-amperes  for  synchronous  generator, 
rated  horsepower  times  746  for  induction  motor 
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ttto,  :  rated  operating  electrical  frequency 


With  these  defined,  impedance  and  tuque  base  quantities  may  be  derived 


(83) 


Use  of  the  pu  system  slightly  changes  the  state  equations  for  die  synchronous 
machine  and  induction  machine.  In  the  pu  system  die  inertia  (7)  is  replaced  by  the  inertia 
constant  (H)  which  has  units  of  seconds.  The  two  inertia  terms  are  related  by 


The  electrical  torque  equations,  (41)  and  (68),  are  altered  in  the  pu  system.  For  the 
synchronous  machine  they  are  written  as 

Tt  -  +  iju  +  ft  - 

—  &)«-*> 

and  for  the  induction  machine  as 


Ti  =  Xyiiqjfr  ~  ijjqr ) 


(86) 
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Use  of  the  pu  system  introduces  a  compatibility  problem  for  the  system  model  Since 
it  is  desired  to  build  a  model  where  sources  and  loads  are  interconnected,  all  submodels 
must  be  referenced  to  the  same  base  values.  The  generator  and  motor  models  share  a 
common  base  voltage  and  electrical  frequency;  however,  base  power  is  different  in  each 
submodel.  For  the  work  presented  here,  the  synchronous  machine  drives  the  selection  of 
base  quantities.  This  requires  that  the  other  submodels  be  referenced  to  the  synchronous 
machine.  By  manipulation  of  (83)  and  (84),  the  machine  parameters  of  the  induction 
motor  may  be  put  in  die  synchronous  generator  base  as  follows: 


7 

pu(sync  _mack) 


^ pm  (sync _  mack) 


r bast  (  sync  _  mack) 

V  ^bast  (  kid  _  mack)  J 
f  p 

r bast ( md  ^  mack) 

^  ^bast(sync  ^mack)  J 


- pu(md_mack ) 


H 


pt(bid  _mack) 


(87) 


C.  SYSTEM  MODEL  IMPLEMENTATION  IN  ACSL 

The  ACSL  program  code  is  contained  in  Appendix  B.  The  description  of  system 
model  implementation  follows  the  layout  of  the  program  code.  In  order  to  make  the 
ACSL  program  easier  to  follow,  throughout  this  section  the  variables  will  be  referred  to 
by  the  same  name  used  in  the  ACSL  code.  The  block  diagram  of  Figure  17  represents  the 
system  simulation  as  implemented  in  ACSL.  The  blocks  labeled  Synchronous  Generator 
Model  (SGM)  and  Load  1  Model  (L1M)  through  Load  n  Model  (LnM)  represent  the 
stand-alone  state-space  models  developed  in  Chapters  II  and  m.  The  models  are 
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Figure  17.  Total  system  rimulition  block  diagram  eg  implemented  In  ACSL. 


implemented  with  current  stales.  At  each  time  step  of  the  simulation,  vectors  representing 
the  model  currents  (j&  and  tf)  and  current  derivatives  (iad  and  ild)  are  computed.  Each  of 
these  models  requires  bus  voltage  as  an  input  The  bus  voltage  to  the  load  may  experience 
a  line  loss. 

The  orthogonal  bus  voltage  values,  vqs  and  vds,  are  generated  by  the  Bus  Voltage 
Model  (BVM).  For  the  system  simulation  presented  here,  the  DASSL  bus  voltage  model 
is  used.  In  order  to  set  up  the  implicit  solution,  the  BVM  requires  all  q  and  d  axis  source 
and  load  current  and  current  derivative  terms  as  inputs.  The  model  forms  a  residual  from 
these  inputs  based  on  KCL.  The  bus  voltage  is  then  computed,  using  the  DASSL  routine, 
at  each  step  of  the  simulation.  The  DASSL  routine  drives  the  residual  close  to  zero. 

The  SGM  requires,  in  addition  to  bus  voltage,  an  input  torque  (71)  and  field  winding 
voltage  (Vfd).  The  SGM  outputs  are  currents  (is)  current  derivatives  (isd)  and  rotor 
speed  (wr).  Speed  control  for  the  SGM  is  provided  by  the  gas  turbine  and  governor  which 
provides  the  input  torque  based  on  the  behavior  of  wr.  The  voltage  control  is 
accomplished  by  the  regulator/exciter  which  uses  the  rectified  magnitude  of  the  Inis 
voltage  to  provide  the  appropriate  level  of  Vfd. 

L1M  through  LnM  receive  bus  voltage  inputs  vql  and  vdl.  These  voltage  quantities 
are  the  BVM  outputs  modified  by  accounting  for  transmission  line  resistance  and 
reactance.  Additionally,  these  load  models  require  an  input  representing  the  electrical 
frequency.  For  this  system  model,  the  electrical  frequency  is  the  nr  output  from  the  SGM. 

The  block  labeled  Inverse  Park's  Transformation  (IPT)  allows  the  system  voltages 
and  currents,  which  are  expressed  in  the  orthogonal  reference  frame,  to  be  changed  to  the 
a-b-c  reference.  The  IPT  uses  thtr,  the  variable  representing  rotor  position,  to  perform 
the  transformation.  This  variable  is  obtained  by  integrating  rotor  speed  with  respect  to 
time. 
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1.  Program  and  Initial  Sections 


The  first  few  hues  of  the  program  set  simulation  parameters.  It  is  here  that 
things  like  integration  algorithm,  integration  step  size  and  communication  interval  for  the 
output  file  may  be  specified.  The  default  integration  algorithm  is  a  Runge-Kutta  fourth 
order  routine.  The  step  size  is  selected,  based  on  the  Nyquist  criteria,  so  that  the  dynamic 
response  transients  can  be  computed  and  printed.  Following  the  setting  of  simulation 
parameters  the  machine  parameters  needed  for  each  state-space  submodel  are  entered. 

The  initial  section  of  the  code  computes  all  coefficients  needed  for  the  source 
and  load  submodels  using  the  machine-  parameters.  Starting  with  the  synchronous 
generator,  the  elements  of  the  inverse  B  matrix  are  computed.  The  expressions  for  these 
elements  are  obtained  from  the  MATHCAD  output  of  Appendix  A.  After  the  inverse 
matrix  elements  are  computed,  the  terms  for  the  linear  and  nonlinear  matrices  are 
calculated.  The  procedure  is  repeated  for  each  induction  motor  model.  Finally  initial 
conditions  for  each  integration  in  the  simulation  are  entered.  These  are  obtained  by  doing 
steady  state  calculations  or  by  running  the  model  with  initial  conditions  set  equal  to  zero 
until  the  model  reaches  steady  state. 

2.  Dynamic  Section 

The  dynamic  section  of  the  ACSL  code  is  the  heart  of  the  simulation.  This 
section  contains  all  the  differential  equations  describing  the  models.  It  is  this  section  of 
the  code  which  is  executed  at  each  time  step  of  the  simulation.  The  speed  of  execution  is 
a  function  of  the  number  of  integrations  which  must  be  performed,  the  integration 
algorithm  selected,  integration  step  size  and  the  total  time  which  must  be  simulated. 
Execution  time  is  also  affected  by  the  DASSL  routine  which  slows  the  simulation  down  at 
each  time  step  by  varying  amounts  in  order  to  drive  the  bus  voltage  to  a  value  that  satisfies 
the  constraint  equation. 
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The  constants  vref  and  wref  are  the  commanded  inputs  for  voltage  magnitude 
and  speed  The  constants  mot_onX  and  cbl  represent  circuit  breakers  for  die  loads.  If 
these  are  set  to  zero,  no  voltage  is  applied  to  the  load  The  induction  motors  each  have  a 
speed  square  law  load  applied  to  them.  This  type  of  load  is  typical  for  a  pump,  a  common 
application  for  induction  motors  on  board  ship.  The  R-L  load  parameters  are  initially  set 
to  large  values  so  that,  unless  they  are  changed  the  load  will  be  very  small  even  with  bus 
voltage  applied  (when  cbl  is  set  to  one). 

Next  some  quantities  are  derived  so  that  they  may  be  output  A  ripple  term  is 
added  to  the  voltage  magnitude  to  simulate  the  output  of  a  rectifier.  To  make  the  output 
more  standard,  torque  quantities  for  the  three  motors  are  changed  from  the  generator 
torque  base  to  the  torque  base  for  each  induction  motor. 

The  exciter  and  prime  mover/govemor  models  come  from  the  s-domain  models 
of  Figures  19  and  20.  The  following  technique  was  used  to  convert  the  transfer  function 
form  of  the  models  to  the  form  seen  in  the  program: 

1  +  v 

B  +  Bxas  -  AKa  (88) 

Bxa  =  AKa  -  B 


where  A  is  the  input,  B  is  the  output  and  B  =  Bs.  Then  this  piece  of  the  larger  system 
may  be  represented  in  ACSL  as 


zl 

B  =  [J  A]  +  h 


(89) 
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This  form,  unlike  the  transfer  function  form,  allows  initial  conditions  to  be  input  for 
integration  if  desired. 

One  of  the  states  for  the  synchronous  machine  is  nr,  the  rotor  speed.  This 
quantity  is  also  the  basis  for  the  reference  frame  chosen  for  die  models.  In  order  to  use 
the  inverse  Park's  transformation  to  convert  q-d-0  quantities  to  a-b-c  quantities,  the  rotor 
position  ( thtr )  must  be  derived.  This  is  done  by  integrating  rotor  speed  with  aspect  to 
time. 

Although  it  has  less  meaning  in  the  case  of  an  isolated  generator  than  in  an 
infinite  bus  or  multiple  generator  study,  rotor  angle  (del)  is  computed.  Figure  18  is  a 
partial  phasor  diagram  which  shows  how  these  quantities  are  related.  The  rotor  angle  is 
sometimes  referred  to  as  the  torque  angle  to  avoid  confusion  with  thtr.  The  physical 
position  of  the  rotor,  thtr,  is  constantly  changing  as  the  prime  mover  drives  the  generator; 
however,  for  a  given  load  on  the  machine  del  is  &  constant 


Figure  18.  Relationship  between  Vigs,  Vds,  Vos  and  del. 

The  state  equation  section  of  th  code  contains  the  ACSL  implementation  of  the 
synchronous  machine  and  induction  machine  equations  in  explicit  form.  The  coefficients 
on  the  right-hand  side  of  the  state  equations  are  named  so  that  they  may  be  easily 
identified.  Coefficients  beginning  with  L  multiply  with  linear  state  terms  Those  with  an  N 
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multiply  with  nonlinear  terms  (those  involving  a  current-speed  product).  The  number 
appended  to  the  name  gives  the  position  of  the  coefficient  in  its  parent  matrix.  In  the 
induction  motor  equations  there  are  two  additional  elements  to  the  naming  conventions 
just  described.  Because  there  are  multiple  motor  models,  a  third  digit  is  added  to  the  end 
of  the  coefficient  name  to  indicate  the  motor  number.  A  distinction  is  also  made  between 
the  nonlinear  coefficients.  Those  involving  the  electrical  speed  only  begin  with  NE,  and 
those  involving  the  slip  speed  or  difference  between  electrical  speed  and  motor  rotor 
speed  begin  with  ND. 

The  DASSL  bus  voltage  model  is  simply  an  expanded  version  of  the  model  used 
in  the  example  of  Chapter  IV.  A  residual  of  the  currents  and  current  derivatives  is  formed 
for  each  q-d-0  axis  .  Then  the  voltage  value  is  obtained  using  the  implicit  system  solver. 
A  line  loss  model  based  on  the  R-L  model  of  Chapter  in  modifies  the  value  of  line  voltage 
applied  to  the  motor  loads. 

The  final  few  lines  of  code  convert  the  voltage  and  current  results  to  the  a-b-c 
reference  frame  for  output  and  set  a  simulation  stop  time.  The  ACSL  system  compiles  the 
model  in  FORTRAN  for  execution.  The  executable  FORTRAN  form  of  the  model  runs 
faster  than  models  produced  in  some  other  simulation  software.  The  main  model  may  be 
linked  to  one  or  more  command  files  which  allow  the  user  to  more  easily  exercise  the 
model  under  a  variety  of  conditions  and  then  obtain  the  desired  output  An  command  file, 
which  was  used  with  this  simulation,  follows  the  system  model  code  in  Appendix  B. 
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VL  SYSTEM  RESPONSE  AND  MODEL  VALIDATION 

The  model  for  the  isolated  power  system  has  now  been  fully  developed.  The  system 
simulation  is  extremely  modular  and  relies  on  accepted  models  for  the  source  and  load 
submodels.  These  accepted  models  for  a  synchronous  generator,  induction  motor  and 
resistive-inductive  load  have  been  extensively  validated. 

Other  submodels  were  needed  to  complete  the  system.  The  DASSL  submodel  for 
the  bus  voltage  was  demonstrated  in  Chapter  IV.  Chapter  V  added  accepted  field  exciter 
and  prime  mover/govemor  models  to  the  system. 

Although  the  separate  pieces  of  the  system  have  been  validated  to  one  degree  or 
another,  the  whole  system  must  be  tested  against  some  standard  for  validation.  Mayer  and 
Wasynczuk  [Ref.  5.]  of  Purdue  University  presented  a  simulation  of  a  portion  of  the  USS 
Arleigh  Burke  (DDG-51)  power  distribution  system  which  will  be  used  as  the  standard  for 
comparison.  This  model  will  be  referred  to  as  the  Purdue  model 

A.  DESCRIPTION  OF  THE  PURDUE  MODEL 

The  Purdue  model  is  a  systematic  method  for  taking  the  differential  algebraic 
equations  describing  a  power  system  and  using  them  to  establish  a  conventional  state- 
space  model.  On  each  bus  of  the  system,  one  machine  is  designated  as  the  root  machine 
and  any  other  connected  machine  is  a  nonroot  machine.  A  root  machine  has  current  and 
current  derivatives  as  its  inputs  and  stator  voltages  as  its  outputs.  The  root  machine 
inputs  are  formed  by  summing  the  currents  and  current  derivatives  from  all  connected 
nonroot  models.  After  establishing  forms  for  the  root  and  nonroot  models,  an 
interconnection  procedure  is  established  based  on  the  KCL  constraint 

The  interconnection  procedure  is  conceptually  similar  to  the  bus  voltage  equation 
development  presented  in  Chapter  IV.  Based  on  the  linearity  of  the  derivative  operator, 

62 


expressions  for  the  current  derivatives  for  the  root  and  nonroot  models  are  set  equal,  thus 
eliminating  the  derivative  terms.  Then,  after  complicated  matrix  algebraic  manipulation, 
an  equation  for  determining  stator  voltages  from  the  states  only  is  produced. 

Mayer  and  Wasynczuk  validate  their  model  by  comparing  it  with  the  output  produced 
by  the  Power  System  Simulator  at  Purdue  University.  This  facility  has  been  used 
extensively  in  system  design  work  and  provides  detailed  three-phase  output  based  on 
state-of-the-art  representations  of  the  power  system  components. 

B.  VALIDATION  BY  COMPARISON  WITH  THE  PURDUE  MODEL 

Mayer  and  Wasynczuk  describe  the  scenario  and  provide  all  the  model  parameters  for 
the  results  presented  in  their  paper.  The  model  parameters  are  contained  in  Table  1.  The 
generator  and  induction  motor  parameters  are  all  per-unit  values  with  a  450V  rms,  line-to- 
line,  base  voltage.  The  base  power  for  the  generator  is  3125  KVA,  and  for  the  motors  is 
determined  by  the  horsepower  rating. 

For  the  comparison  simulation,  the  generator  is  initially  in  a  steady  state  unloaded 
condition.  It  is  under  closed  loop  regulation,  operating  at  rated  voltage  and  speed  with 
stator  currents  zero.  At  an  arbitrary  time,  a  circuit  breaker  is  closed  energizing  all  three 
induction  motors.  The  three  motors  draw  large  start-up  currents  which  cause  bus  voltage 
to  dip  initially  before  the  voltage  regulator  and  exciter  circuit  can  react  and  return  the  bus 
voltage  to  the  commanded  magnitude.  The  initial  large  currents  also  produce  a  large 
electrical  torque  in  the  generator  which  tends  to  slow  rotor  speed.  The  prime 
mover/govemor  reacts  to  the  speed  change  by  applying  more  input  torque  to  return  the 
system  to  commanded  speed.  Most  of  the  system  transient  behavior  is  complete  in  three 
seconds.  Figures  19  and  20  compare  the  results  of  the  DASSL  based  simulation  with  the 
results  obtained  by  the  Purdue  model. 
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TABLE  1.  SYSTEM  MODEL  PARAMETERS 


Prime  Mover/Governor 

It.  =  22.5 

T„  =  0.55 

■VI 

wmm 

_ Wp,n.  =  0.23 

ll  _ 

ter 

msmm 

i 

To,  =0.15 
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B  =  0.3 

I  Synchronous  Generator 
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hhhh 

h^ssiis 
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Induction  Motors 
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Hp 

200 

150 

40 
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x_ 
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X. 

0.0655 
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0.0587 

IV 

0.0261 

0.0165 

0.0165 

H 

0.922 

1.524 

1.054 
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Figure  19.  Comparison  of  results:  terminal  voltage  (vj),  field  voltage  (Ex),  prime 
mover  torque  (Tpnih  ™tor  speed  (<%■),  phase  A  voltage  (V^;  (a)  DASSL  model, 

(b)Purdue  model. 
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Figure  20.  ConqNuisQDof  remits:  phase  A  current  (ig),  motor  1  torque  (Tej),  motor 
1  speed  ( (Of  j), motor  2  torque  (Te2)>  motor  2  speed  (tori)*  (®)  DASSL  model, 

(b)Purdue  model. 
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The  DASSL  based  model  is  in  excellent  agreement  with  the  Purdue  model  results. 
Additionally  it  agrees  well  with  the  expected  behavior  of  such  a  system.  Bus  voltage 
transient  behavior  is  tracked  (hiring  dynamic  loading  of  the  system.  During  the  simulation 
the  difference  between  generator  current  and  the  sum  of  the  load  currents  stayed  bounded 
in  the  micro  amp  range. 

Both  simulations  were  implemented  in  ACSL  using  a  fourth-order,  Runge-Kutta 
integration  algorithm.  Using  an  integration  step  size  of  1.0  millisecond,  the  DASSL  based 
model  required  5.4  seconds  of  cpu  time  on  a  Sun  SPARC  10  Model  41  workstation.  For 
comparison,  the  generator  and  three  induction  motor  simulation  using  an  infinite  bus 
voltage  model  used  2.1  seconds  of  cpu  time  on  die  same  workstation. 

C.  DASSL  MODEL  RESPONSE  WITH  UNBALANCED  LOAD 

In  three-phase  power  systems  every  effort  is  made  by  the  system  designers  to  keep 
the  phase  loads  equal.  In  practice  this  is  impossible.  All  three  phase  systems  experience 
some  degree  of  unbalanced  loading.  On  board  a  ship  this  is  a  common  problem, 
particularly  on  older  ships.  Partial  grounds,  lighting  alterations  and  equipment  updates 
are  some  factors  that  contribute  to  the  problem.  Behavior  of  the  system  presented  in  the 
previous  section  to  an  unbalanced  loading  condition  will  be  investigated. 

An  unbalanced  R-L  load  in  parallel  with  the  three  induction  machines  simulates  a 
situation  which  frequently  occurs  on  board  ship.  Single  phase  lighting  loads  are  served  by 
tapping  a  lighting  circuit  transformer  primary  into  one  phase  of  the  three  phase  power 
system.  This  often  results  in  different  loading  conditions  on  the  three  phases. 
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I.  Hie  Unbalanced  Load 

Unbalanced  loading  introduces  significantly  nunc  complexity  to  the  state 
equations  in  the  q-d-0  reference  frame.  In  the  R-L  load  development  of  Chapter  HI, 
hflianrad  loading  was  assumed.  If  the  phase  resistance  values  are  allowed  to  be  unequal 
the  resistance  matrix  becomes 


*i  = 


r«i  o 

0  rw 
0  0 


0 

0 

ra. 


The  transformation  to  the  orthogonal  reference  frame  results  in  the  Ktr,(K,~l)  matrix 
taking  the  form 


cos2  A  +  rbl  cos2  B 
■•■+rbl  cos2  C 

cos  A  sin  A  +  rbl  cos  B  sin  B 
•  ••  +rbl  cos  C  sin  C 

■j  (r^  cos  A  +  rN  cos  B 
•••  +rbl  cos  C) 


cos  A  sin  A  +  rwcos£sin/J  cos  A  +  rw  cos  B 


—  +ij|  cos  C  sin  C 

^  sin2  A  +  sin2  fl 
•••  +ru  sin2  C 

^  sin  A  +  rbl  sin  B 
•  •+rwsinC) 


•••  +/^  cos  C 

sin  A  +  rw  sin  B 
•••  sin  C 


^2  (r«i  +  ru  +  'ct) 


(90) 


where 


a  =  er 
‘•••-T 

c.,.f 
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This  significantly  complicates  the  stale  equations.  All  three  states  now  appear  in 
each  of  the  three  equations  describing  the  R-L  load  in  the  q-d-0  reference  frame.  Note 
that  an  unbalanced  inductance  matrix  will  cause  all  three  state  derivative  terms  to  appear 
in  all  state  equations  (an  implicit  form).  Putting  such  a  system  of  equations  in  explicit 
form  will  not  be  dealt  with  here. 

2.  Simulation  Results  with  Unbalanced  Loading 

For  the  unbalanced  load  study,  the  model  is  started  and  loaded  in  an  identical 
manner  to  the  simulation  of  Figures  19  and  20.  Once  the  system  is  in  steady  state,  a 
circuit  breaker  is  closed  which  connects  an  unbalanced  R-L  load  in  parallel  with  the  motor 
loads.  The  per-unit  unbalanced  load  parameters  used  are 

-  5.0 
ru  =30.0 
rel  =  5.0 
Xt  =  3.0 

After  application  of  the  unbalanced  load,  the  system  is  again  allowed  to  come  to 
steady  state.  The  effect  of  this  load  on  the  phase  currents  may  be  seen  in  Figure  21.  The 
phase  currents  are  visibly  unequal  The  magnitude  of  is  about  half  the  magnitude  of  the 
other  phase  currents.  Of  more  interest  is  the  manner  in  which  other  variables  are  affected. 

Figure  22  is  a  plot  of  several  variables  affected  by  the  unbalanced  load.  The 
unbalanced  load  causes  oscillations  to  occur  in  the  generator  electrical  torque.  These 
variations  in  turn  cause  the  rotor  to  oscillate.  The  field  winding  current  also  exhibits  this 
oscillatory  behavior  due  to  the  motion  of  the  rotor.  Not  only  is  the  generator  affected,  but 
the  plot  of  torque  produced  by  induction  motor  two  shows  oscillations.  Such  oscillations 
can  be  potentially  damaging  to  equipment 
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Figure  22.  System  response  to  unbalanced  load;  generator  electrical  torque  (7^), 
generator  rotor  speed  (av/Q9&),  field  current  (fa)  and  motor  2  torque  (Tl2). 
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D.  DASSL  BASED  MODEL  FLEXIBILITY 

Many  other  types  of  studies  may  be  done  with  the  DASSL  based  model  The  system 
is  extremely  flexible.  For  example,  the  load  on  one  or  more  of  die  induction  motors  .  iuld 
be  suddenly  changed.  The  full  transient  effect  of  this  change  could  then  be  observed 
throughout  the  system.  The  response  of  the  generator,  generator  control  systems  and 
other  loads  could  all  be  studied.  This  type  of  capability  is  extremely  valuable  to  the 
designers  of  isolated  power  systems. 
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m  CONCLUSIONS  AND  FUTURE  WORK 


There  is  a  need  for  power  system  simulations  which  do  not  rely  on  the  infinite  bus 
voltage  assumption.  In  die  recent  post,  the  limited  capabilities  of  computer  hardware  and 
software  mack;  tire  modeling  of  isolated  power  systems  difficult  Designers  relied  on 
questionable  equations  to  approximate  terminal  voltage  dip  under  various  loading 
conditions.  Reduced  order  modeling  was  done  which  provided  some  data  at  the  expense 
of  losing  transient  behavior  results.  It  has  been  demonstrated  that  by  treating  the 
equations  for  tire  power  system  as  a  set  of  differential  algebraic  equations  and  using  a 
proven  DAE  solver,  excellent  results  can  be  achieved. 

A.  ADVANTAGES  OF  THE  DASSL  MODEL 

The  approach  presented  in  this  diesis  has  several  advantages  over  methods  that  have 
been  used  in  the  past.  Some  of  the  advantages  are: 

«  the  system  model  is  highly  modular  in  design  (submodels) 

•  the  DASSL  bus  voltage  submodel  constraint  equation  makes  the  model  simple 
to  expand 

•  the  model  provides  transient  data 

•  the  submodels  are  standard,  well  validated  state-space  models 

•  simulation  speed  is  excellent 

«  the  model  uses  a  highly  regarded,  commercially  available  DAE  solving  routine 

B.  DISADVANTAGES 

There  are  some  problems  with  the  system  model  which  still  must  be  overcome. 
Some  of  the  possible  disadvantages  and  limitation  of  the  model  are: 
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•  the  source  and  load  submodels  need  to  be  developed  with  current  states  and  in 
the  q-d-0  reference  frame  (not  always  the  most  convenient  form) 

•  die  total  system  model  has  not  been  validated  against  hardware  results 

•  the  generator  and  motor  models  do  not  include  saturation  effects 

•  the  generator  and  motor  models  assume  sinusoidally  distributed  windings 

The  last  disadvantage  is  ready  a  disadvantage  of  the  standard  development  of  the 
synchronous  and  induction  machine  models.  The  fact  that  the  windings  are  actually  not 
perfectly  distributed  introduces  harmonics  on  the  system  bus.  Generally  these  effects  are 
minor  and  may  be  neglected;  however,  some  types  of  loads  (solid  state  power  converters 
for  example)  may  be  highly  sensitive  to  these  neglected  harmonics. 

C.  FUTURE  WORK 

In  order  to  realize  a  benefit  from  the  DASSL  model  approach  significant  work 
remains  to  be  done.  One  of  the  most  important  tasks  that  could  be  accomplished  is  a 
hardware  validation  of  the  model.  Figure  23  is  a  block  diagram  of  a  possible  hardware 
configuration  for  accomplishing  validation. 


Figure  23.  Possible  hardware  configuration  for  model  validation. 


This  configuration  has  the  advantage  of  being  complex  enough  to  validate  the 
DASSL  model  and  simple  enough  to  build  and  test  in  the  lab.  Some  other  suggestions  for 
future  work  with  this  basic  model  are: 

•  addition  of  one  or  more  parallel  generators 

•  development  of  other  types  of  load  models 

•  include  winding  saturation  effects 

•  try  other  DAE  solution  methods  which  offer  the  promise  of  improved 
efficiency  (for  example,  Halin  [Ref.  12]  reports  a  significant  improvement  over 
DASSL  with  his  method) 

This  list  is  by  no  means  conclusive  and  much  work  remains  to  be  done  in  this  area. 

The  process  of  engineering  design  involves  trade-offs.  The  DASSL  model,  because  it 
involves  multiple  iterations  at  each  simulation  time  step,  requires  more  than  double  the  cpu 
time  than  that  used  by  the  infinite  bus  model.  At  one  time  this  cost  may  have  been 
considered  too  great  for  the  benefit  derived.  However,  today's  simulation  capabilities 
make  the  DASSL  based  finite  bus  model  a  practical  design  tool,  worthy  of  continued 
investigation. 
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APPENDIX  A:  CONVERTING  STATE  EQUATIONS  TO  EXPLICIT  FORM 


In  Chapters  Q  and  m  the  state  equations  for  the  synchronous  machine  and  the 
induction  machine  were  developed  as  implicit  equations  with  the  currents  as  states.  Those 
equations  had  the  form 

v  =  ALi  +  Ano )/  +  B pi 

where  Al  is  the  linear  terms  matrix  and  AN  is  the  nonlinear  terms  matrix. 

In  order  to  convert  the  equations  to  implicit  form  the  B  matrix  must  be  inverted  so 
that  the  state  derivative  vector  may  be  isolated  on  the  left-hand  side  of  the  equation.  Since 
numerical  matrix  inversion  can  often  lead  to  problems  in  the  case  of  poorly  conditioned 
matrices,  the  matrix  inversion  for  the  two  models  was  done  using  the  MAPLE  symbolic 
engine  in  MATHCAD  4.0. 

THE  SYNCHRONOUS  MACHINE 


-b  0  0  •  0  0 

0  -k  0  0  b  b 

-1 

where  the  following  definitions  apply: 

0  0  -c  0  0  0 

-i  0  0  g  0  0 

8  =  b  =  ;  c  =  X„ 

b2  e  b2 

o  -Z-  0  0  -b  — 
d  d  d 

0  -b  0  0  b  f 

d  =  RM;e  =  XM;f  =  XM 
g  =  XM;b  *  Xq;k  =  X„ 
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THE  INDUCTION  MACHINE 


B-' 


X„  0  0  X,,  0  0 

0  X.  0  0  X||  0 

0  0  \  0  0  0 

Xj,  0  0  xn  0  0 

0  X,,  0  0  X„  0 

0  0  0  0  o  xk 
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THE  STATE  EQUATIONS  FORMED  EXPLICITLY  FOR  ACSL 

The  final  form  of  the  state  equations  isolates  the  current  derivative  vector  on  the  left- 
hand  side.  In  the  ACSL  program  code  an  equation  is  written  for  each  current  derivative. 
In  order  to  develop  an  equation  fc-  each  state  derivative,  the  matrix  equations  must  be 
multiplied  out.  This  development  will  be  demonstrated  for  the  synchronous  machine,  the 
induction  machine  is  treated  in  an  identical  manner.  Equation  (42)describes  the  procedure 
in  compact  from. 


pi  =  Li+N(Dr£+ Vy 
V  =  B_I 
L  =  V(-At) 

N  =  V(-An) 


Using  MATHCAD  again,  the  state  equations  may  be  developed  as  seen  in  the  ACSL 
code  of  Appendix  B.  First  the  L  and  N  matrices  must  be  computed 
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0 

r 

0 

0 

0 

0 

0 

V22 

0 

0 

V25 

V26 

• 

0 

0 

V33 

0 

0 

0 

0 

0 

r. 

0 

0 

0 

V4l 

0 

0 

V44 

0 

0 

0 

0 

0 

0 

0 

0 

V52 

0 

0 

V35 

VS6 

0 

0 

0 

0 

x-d 

0 

0 

V62 

0 

0 

VOS 

V66 

0 

0 

0 

0 

0 

-rM 

L  = 


VUr,  0  0  -VMt^  0  0 

0  V22r,  0  0  -V25-X—  -V26 rM 

0  0  V33rt  0  0  0 

V41rt  0  0  -V44-fti|  0  0 

0  V52r§  0  0  -V55  X_  -V56-ru 

0  V62r,  0  0  -VMX^  -V«-rM 
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VII  0  0  VI4  0  0 

0  V22  0  0  V23  V2» 

0  0  V33  0  0  0 

V4I  0  0  V44  0  0 

0  VJ2  0  0  V55  VJ4 

0  V62  0  0  V«5  V«6 


0—00 
A  0  o5s 

•k  % 

0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 


„  ** 

X* 

0 

VII — 

0 

0 

-vu  — 

-YU¬ 

•S 

•k 

'S 

X, 

-V22 — * 

0 

0 

V22^ 

0 

0 

®k 

■k 

0 

0 

0 

0 

0 

0 

x< 

X* 

0 

V41 — 

0 

0 

-V4I  — 

-V41— — 

®k 

“k 

®k 

-VS2-1 

0 

0 

X„ 

V32-21 

0 

0 

®k 

“k 

-V62—-S- 

0 

0 

V02-— t 

0 

0 

®k 

®k 

then  the  terms  of  the  equation  are  multiplied  out 


Lit  0  0  U4  0  0 

0  U2  0  0  L23  L26 

0  0  L33  0  0  0 

Ml  0  0  M4  0  0 

0  LSI  0  0  L55  L56 

0  M2  0  0  US  M6 


L22'iA  «-L25-^  +  L2fr^ 
L33^, 

Ml^+M4^ 
LS2\  i-LSS  hrLS&\d 
L62  ^  ■►MS  ^ +  L66-iu 


No),i  = 


0  N12  0  0  N1S  NI6  ® 

N21  0  0  N24  0  0  '* 

0  0  0  0  0  0  S 

0  N42  0  0  N4J  N46  *  ^ 

NSI  0  0  NS4  0  0  ^ 

N61  0  0  N64  0  0 


N12  ®, *  NIS®,^  N16  0»r  ^, 

N2l«r^  +N24«r-^ 

0 

N42®f *  l*45®,'  W  ♦  N46'®, 
N51e»(i^  +-N54®t^ 
N61®r  i^  ♦N04  «f\( 
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V*- 


fxnally,  the  expression  for  the  state  derivative  vector  may  be  written  as  the  sum  of  the 
three  terms 


NI2  ♦  NlS  fflj  ij,  +  N16  «r +  Lit ♦  L14 ♦  Vllv^ 

N21  +  N24oBfibj  +  L22ifc  + L^S-^  +■  L26-^  ♦  V22vfc  +  V25vM 

'** 

N42«0fi(fc  +  N45a»r +  N46a>f +  Uli^  +  LM-i^  +  V41v^ 

N51  ayi  ♦  N54a>rib}  +  L52 +■  LSI ^  + 1^6^  +  V52-V,,  +  V55vM 

*u 

N61  <Dfi  +  -t- 1^2-^  +  L65 +  L66 +•  +  V65vfd 
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APPENDIX  B:  ACSL  CODE 


A.  BUS  VOLTAGE  EQUATION  MODEL 

!l!!!!l!!IU!lll!t(lllll(l!III!l!!l!tfif!!!li{!IIII!!!!!l!tl!lltli!!ill!!llllll!lll!llll!l!llll!llllllltlllltill 

II 

II  Mark  Kipps  NPS  Monterey 

II 

II  Program  to  demonstrate  validity  of  bus  voltage  equrtion  modal  . 
II  Example  circuit  is  solved  using  two  methods  and  the  results 
II  are  compared. 

II 

l!lll!llll!!llllllllllflllll!llllll!lltlllltltllllll!lllil!!lllllll!lllllllt!!llt!lllllllllltllllll!ll!l!llllll! 

PROGRAM 

NSTEPS  nstp*1 
CINTERVAL  tint .  1e-2 
MAXTERVAL  maxt-  1e-3 

DYNAMIC 

DERIVATIVE 

1-CircuX  parameters 
CONSTANT  R1-1.0 
CONSTANT  R2<«5.0 
CONSTANT  LI  -0.6 
CONSTANT  L2-0.2 

(--Source  voltage 

Vs  -PULSE(0.0, 2.5, 1.25) 

1-State  equations  for  one  loop  solution 

id  .  -(R1  +R2)/(L1  +L2)*i  +  Vs/(L1  +L2) 

1  -  INTEG(kJ,0.0) 

Vnode  -  L2*id  +  R2*i 

1-State  equations  for  the  two  sub-model  solution 

lid  -  -R1/L1*I1  +  Vs/Ll  -  WL1 

il  -  INTEG(i1d,0.0) 

i2d  « -R2/L2*i2  +  WL2 

2  m  INTEG(i2d,0.0) 
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I-Bui  voftage  equation  for  aub-modais 

VI  -fl.1*L2)*(-R1H/L1  ♦  R212/L2  +  V*Liy(L1  +  L2) 

l-Dfferencee  tor  ouq>ut 

detv  *  Vt  -  vnode 
del  -11-12 
defid  -ild-Bd 

END  1  of  derivative 

CONSTANT  tetop  -  20. 

TERMT(t  .QE.  tetop) 

END  !  of  dynamic 
END  I  of  program 
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B.  DASSL  BUS  VOLTAGE  MODEL 


!!t!lll!lll!ttttlItt!!!t!tfl!f!tltlllt!!l!imil!IIII!!l!l!f!!!!!!!l!lt!tlll!!IIIt!!l!l!!tlll!!!!ltlttllll!tt!t! 

I! 

!!  Marie  Kjppe  NPS  Monterey 

!! 

!!  Program  to  demonstrate  vattdty  of  DASSL  bus  voltage  model.  !! 
II  Example  circuit  is  solved  using  two  methods  and  the  results 
II  are  compared. 

II 

IfllKlilllllflllilllflflllllllllillllllltlflllllillllHIHIIIIIIIfltllllHflltlllflllllllHIIIIIIIIIIIIIIHItll 

PROGRAM 

NSTEPS  nstp  .  1 
CINTERVAL  dnt  -  1e-2 
MAXTERVAL  maxt  »  1e-3 

DYNAMIC 

DERIVATIVE 

! --Circuit  parameters 
CONSTANT  R1-1.0 
CONSTANT  R2-5.0 
CONSTANT  LI-0.6 
CONSTANT  L2-0  2 

(-Source  voltage 

Vs  =PULSE(0.0, 2.5, 1.25) 

! --State  equations  for  one  loop  solution 

id  -  -(R1  +R2y(L1  +L2)*i  +  Vs/(L1  +L2) 

1  =  INTEG(id,0.0) 

Vnode  »  L2*id  +  R2*i 

!--State  equations  for  the  two  sub-model  solution 

•Id  -  -R1/LH1  +  Vs/Ll  -  WL1 

il  -  INTEG(i1d,0.0) 

i2d  -  -R2/L2*i2  +  WL2 

2  .  INTEG(i2d,0.0) 

!-- DASSL  bus  voltage  model  sums  current  at  the  node 

resi  » lid  + 11  -  i2d  - 12 

Vt  -  IMPLC(resi,0.0) 
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(-DVforancee  tor  oufcx* 

deN  ■  Vl -vnode 
del  —11-12 
deid  -  lid  -  Bd 

END  !  of  derivative 

CONSTANT  tetop-20. 
TERMT(t  QE.  tatop) 


END  I  of  dynamic 
END  I  of  program 


C.  BUS  VOLTAGE  EQUATION  MODEL  UNDER  CONTROL 

!il(llf!tllit!1!lllll!!!l!l!ll!lll!ll!ll!l!l!llli!!!!llfllllll!l!!!!!!t!lttlt!!!l!!l!!lll!lllllllllllllll!!ll!!!  ! 
II  I 

II  Mark  Kipps  NPS  Monterey  I 

II  I 

II  Program  to  demonstrate  validity  of  bus  voltage  equation  model.  II 
I!  Example  circuit  Is  solved  using  two  methods  and  the  resuite  l 

II  are  compared.  System  under  cascade  voftage  control.  I 

II  I 

lll!lll!ll!lllll!lll!lll!!!l!lllllllillll!lllllf!!i!ll!illlillllllll(ltlllllll(llll!lflflllt(flll!fflfllfllll!fl  I 

PROGRAM 

NSTEPS  nstp  =  1 
CINTERVAL  dnt«5e-3 
MAXTERVAL  maxt  =  1e-3 

DYNAMIC 

DERIVATIVE 

l~Circuit  parameters 
CONSTANT  R  1-1.0 

CONSTANT  R2-5.0 
CONSTANT  LI  =0.6 

CONSTANT  L2-0.2 

Vref  =STEP(0.2) 

1-Cascade  voltage  controller  from  root-locus  design 

verr  =Vref-Vt 

void  =200*verr  -  30*v01 

vOI  =INTEG(v01d,0.0) 

Vsd  =v01d  +  10*v01  -  .01  *Vs 

Vs  =INTEG(Vsd,0.0) 


!-- State  equations  for  two  sub-model  solution 

ild  =  -R1/L1*i1  +  Vs/Ll  -  WL1 

11  =  INTEG(lld.O.O) 

i2d  =  -R2/L2*i2  +  WL2 

12  » INTEG(i2d,0.0) 

1-Bus  voltage  equation 
resi  =  ild  +  il  -  i2d  *  i2 

Vt  =(L1  *L2)*(-R1  *i1/L1  +  R2*i2/L2  +  Vsfl.1  )/(L1  +  L2) 


86 


!-Tran«fer  function  form  of  the  system 


ve  »Vref  •  Vb 

DIMENSION  p(3),  q(4) 

CONSTANT  p- 1 .0. 35.0,  250.0,  q- 1 .0, 37.51 , 225.375, 2.25 
Vb  -TRAN(2,3,p,q,50.0*ve) 

(-Differences  for  output 

del  « il  -  i2 
deNd  •  ild  -  i2d 
detv  =Vt- Vb 

END  I  of  derivative 

CONSTANT  tstop  >1.0 
TERMT(t  ,GE.  tstop) 

END  I  of  dynamic 
END  I  of  program 
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D.  DASSL  BUS  VOLTAGE  MODEL  UNDER  CONTROL 

!!flllll!!!!l!f!!!flft!llllf!l)!t!!lltl!fl!lll!!lll!l!ll!ll!ltt(t!t!t!!lttttllltll!lll!(llll!lllllltl11!",a.!1ll  I! 
I!  II 

II  Mark  Kippe  NPS  Monterey  II 

II  II 

I!  Program  to  demonstrate  vaMdty  of  DASSL  bus  voltage  model.  II 

II  Example  circuit  is  solved  using  two  methods  and  the  results  II 

I!  are  compared.  System  under  cascade  voltage  control.  II 

I!  I! 

Illl!ll!!!ltlll!lillli!l!!llllil!llll!lllll!llilllllllllll!!lll!llllll!lllllllllll!iitlill!!!!lli!lllll!tlllllll  H 

PROGRAM 

NSTEPS  nstp  *  1 
CINTERVAL  cint «  5e-3 
MAXTERVAL  maxt=1e-3 

DYNAMIC 

DERIVATIVE 

I— Circuit  parameters 
CONSTANT  R1=1.0 
CONSTANT  R2=5.0 
CONSTANT  LI  =0.6 

CONSTANT  L2=0.2 


Vref  =  STEP(O^) 

!-Cascade  voltage  controller  from  root-locus  design 

verr  =Vref-Vt 

void  =200*verr  -  30*v01 

vOI  =INTEG(v01  d,0.0) 

Vsd  =v01d  +  10*v01  -  .01*Vs 

Vs  =INTEG(Vsd,0.0) 


! --State  equations  for  two  sub-model  solution 

ild  =  -R1/L1*i1  +  Vs/Ll  -  Vt/LI 

11  =  INTEGKild.O.O) 

i2d  =  -R2/L2*i2  +  WL2 

12  =  INTEG(i2d,0.0) 

! --DASSL  bus  voltage  model 
resi  =  ild  +  il  -  i2d  -  i2 
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Vt  .  IMPLC(reei.O.O) 

{-Transfer  function  form  of  the  system  for  comparison 

ve  »Vrsf  -  Vb 

DIMENSION  p(3),  q(4) 

CONSTANT  p-  1 .0, 35.0,  250.0,  1 .0. 37.51 , 225.375, 2.25 

Vb  -TRAN(2.3,p,q,50.0*ve) 

(-Differences  for  output 

del  a  II  *  12 

delid  =  ild  -  i2d 

detv  =  Vt- Vb 

END  !  of  derivative 

CONSTANT  tstopal.O 
TERMT(t  .GE.  tstop) 

END  I  of  dynamic 
END  (of program 
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E.  TOTAL  SYSTEM  MODEL 


Ulff!!lflfll!HI!!fl!HI!lllllfH!!l!ll!!!HllllfH!lllll!lllll!fllfl!!l!!!!!!!!HHH!!l!lf!!IHII!H!H!ll!!l  II 

II  « 

II  MarkKIppe  NPS  Monterey  II 

II  » 

II  Program  to  simulate  a  synchronous  generator  and  loads  on  a  flrtfe  II 

II  system  bus.  Bus  voltage  is  computed  using  the  impicit  equation  II 

I!  solving  routine  OASSL.  II 

I!  H 

l!l!l!!!IIll!!llll!lll!!!!lit!!i!!illliltlll!llfllfHHIIIIHIIillll!!illl»l»IHIIII!ltllli!IIHI>llll>IDIill  H 

PROGRAM 

NSTEPS  nstp  *  1 

CINTERVAL  dnt » 1e-3 

MAXTERVAL  maxl  =  1e-3 


INITIAL 


pi  =  4.0*atan(1.0) 

wb  *  120.*pi 

!!!!!!!!! — Synchronous  Machine  Parameters 


H 

=  2.137 

smpb 

=  3125. 

zb 

=  1.0 

Xs 

*.08'zb 

Xmq 

=  1.0*zb 

Xmd 

=  1.768*zb 

XIKd 

=  .33383*zb 

XNd 

=  13683*zb 

Xlkq 

=  .3298*zb 

Rfd 

=  .00111*zb 

Xkd 

=  Xlkd+Xmd 

Xfd 

=  Xlfd+Xmd 

Xkq 

=  Xikq+Xmq 

Rs 

=  .0051 5*zb 

Rkd 

=  .02397*zb 

Rkq 

=  .061 3*  zb 

Xd 

=  Xs+Xmd 

Xq 

=  Xs+Xmq 

!!!!!!!!! — Induction  Motor  Parameters  200 hp  machine 
impbl  *  149.14 
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zbjml  -smphflmpbl 
Rsjml «  0.01*zbjm1 
Xtojml  •  0.0655*zbjml 
XmJml  » 3.225*zbjml 
Xlr_im1 «  0.0655*zbjm1 
Rr_im1  ■  0.0261 *zbjm1 
HJml  -0.902*140.14/3125. 

Xasjml  -XJsJml+XmJml 
Xrrjml  «  Xlr_im1  +Xm_lm1 

1!!II!!!I — Induction  Motor  Parameters  150hp  machine- 

impb2  -111.85 
zb Jm2  «  smpb/impb2 
Rs_hn2  -  0.0051  *zbjm2 
Xtejm2  »  0.0S53*zbJm2 
Xm_im2  »  2.678*zb_im2 
Xlr_im2  -  0.0553*zbjm2 
Rr_im2  «  0.0165*zbjm2 
H_im2  -1.524*111.85/3125. 

Xss_im2  -  XlsJm2+XmJm2 
Xrr_im2  «  Xlr_im2+Xm_im2 

!!!!!!!!!— Induction  Motor  Parameters  40hp  machine- 

impb3  -  29.83 
zb_im3  -  smpb/impb3 
Rsjm3«  0.005*zb_im3 
Xte_im3  -  0.0587*zb_im3 
Xm_im3  -  2.952*zbjm3 
Xlr_im3  -  0.0587*zb_im3 
Rrjm3  •*  0.0165*zb_im3 
H_im3  - 1.054*29.83/3125. 

Xss_im3  =  Xte_im3+Xm_im3 
Xrr_im3  -  Xlr_im3+Xm_im3 

!!!!!!!!!— Field  Exciter  Parameters - 


Kaa-400. 

Tae  -  0.01 
Kfe  -  0.01 
Tfel  -0.15 
Tfe2  -  0.06 
Kee- 1.0 
Tee  -  0.1 

!!!!!!!!! — Prime  Mover  and  Speed  Governor  Parameters 


Kc  -  22.5 
Tc  -  0.55 
Tfv  -  0.01 
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TO -0.05 
WflOs  -  0.23 
C2gt  -  0.251 
Clgt- 1.3623 
Cgngt-0.5 


llltlllll — Synchronous  Machine  modal  coefficient* - 

! — Inverse  matrix  (obtained  from  MATHCAD) 

V1 1  -  wt>*Xkq/(Xmq**2  -  Xq'Xhq) 

V14  -  -\M>*Xmq/(Xmq**2  *  Xq'Xhq) 

V22  -  wb*(Xmd**2  -  Xfd*XtaJy  & 

(XcTXWXkd  +  (2*Xmd  •  Xkd  -  Xd  -  Xfd)*XmcT2) 

V2S  -  wb*Rfd*(Xkd  •  Xmdy  & 

(Xd*Xfd*XKd  +  (2*Xmd  -  Xkd  •  Xd  -  Xfd)*XmtT2) 

V26  -  wb*Xmd*(Xfd  -  Xmdy  & 

(X<rXM*Xkd  +  (2*Xmd  -  Xkd  -  Xd  -  Xfd)*XmcT2) 

V33  -  -wtVXs 

V41  «  wb*Xmq/(Xmq**2  -  Xq*Xkq) 

V44  -  -wb*Xq/(Xmq**2  -  Xq*Xkq) 

V52  -  wb*Xmd*(Xmd  -  Xkd)/  & 

(XcTXfd*Xkd  +  (2*Xmd  -  Xkd  -  Xd  -  Xfd)*Xmd~2) 

V55  -  wb*Rfd*(Xd*Xkd  -  Xmd"2y  & 

(X<rXfd*Xkd  +  <2*Xmd  -  Xkd  -  Xd  -  XM)*Xmd**2yXmd 
V56  *  wb*Xmd*(Xmd  *  Xdy  & 

(Xd*Xfd*Xkd  +  (2*Xmd  -  Xkd  -  Xd  -  XM)*Xmd“2) 

V62  -  wb*Xmd*(Xmd  -  Xfdy  & 

(Xd*Xfd*Xkri  +  (2*Xmd  -  Xkd  -  Xd  -  X!d)*Xmd"2) 

V65  -  wt>*Rfd*(Xmd  -  Xdy  & 

(Xd*Xfd*Xkd  +  (2*Xmd  -  Xkd  *  Xd  -  Xfd)*Xmd**2) 

V66  -  wb*(XM*Xd  -  Xmd"2y  & 

<Xd*Xfd*Xkd  +  (2*Xmd  -  Xkd  *  Xd  -  Xfd)*Xmd**2) 

!— Linear  terms  matrix 

L11  -  VII'R. 

L14  » -V14*Rkq 
L22  »  V22*Rs 
L25--V25*Xmd 
L26  -  -V26*Rkd 
L33- V33*Rs 
L41  -  V41*Rs 
L44  « -V44*Rkq 
L52-V52*Re 
L55  -  -V55*Xmd 
L56  -  -V56*Rkd 
L62  -  V62*Rs 
L65  -  -V65*Xmd 
L66  -  -V66*Rkd 

(—Nonlinear  terms  matrix 
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N12-Virxdftrt> 

N15  w-VirXmdftib 
N16a-virXmc¥*b 
N21  -  -V22*Xqtob 
N24  -  V22*Xmq/wb 
N42«V41*XdArt> 

N45  >  -V41*Xmd/wb 
N46  -  -V41*Xmd/wb 
N51  -  -V52*Xq/V4) 

N54  «  V52*Xmq/wb 
N61  -  -V62*Xq/wb 
N€4  »  V62*Xmq/wb 

!!!!!!!!! — Coefficients  for  Induction  Motor  Loads- 


MACRO  imcoef  (x) 

I— Inverse  matrix 

Bll&x  =  XnJm&x/(X8sJm&x#XrrJm&x  -  Xm_imAx**2)*wb 
B14&X  *  -Xm_im&x/(Xss_im&x*Xrr_im&x  -  Xm_im&x~2)*wb 
B22&X  -  Xrr_im&x/(X8S_inr>Ax*Xrr_im&x  -  Xm_im&x**2)*wb 
B25&X  « -Xm_im&x/(X88_im&x*Xrr_im&x  -  Xmjm&x"2)*wb 
B33&X  *  1/Xl8jm&x*wb 

B41&X  =  -XmJm&x/(XssJm&x*Xrr_im&x  -  XnUm&x**2)*wb 
B44&X  *  Xss_im4x/(XssJm&x*Xrr_lm&x  -  Xm Jm&X**2)*wb 
B52&X  » -Xm_imAx/(XssJmAx*XrT_im&x  -  Xm_imta~2)*wb 
B55&X  =  Xss_im&x/(X88jm&x*Xrr_im&x  -  Xm _imta**2)*wb 
B66&X  «  1/Xlr_im&x*wb 

LM1  l&x  *  -B11&x*R8_im&x 
LM14&X  «  -B14&x*Rr_im&x 
LM22&X  =  -B22&x*Rs_imAx 
LM25&X  *  -B25&x*Rr_im&x 
LM33&X  =  -B33&x*Rs_im&x 
LM41&X  *  -B41&x*RsJm&x 
LM44&X  -  -B444x*RrJm&x 
LMS2&X  -  -B52&x*Rs_im&x 
LM55&X  -  -B55&x*Rr_im4x 
LM66&X  - -B66&x*RrJm*x 


NE12&X  a  -B1 1  &x*X88_im&x/wb 
NE1 5&x  >  -B1 1  &x*Xm_im&xftvb 
NE21&X  >  B22&x*Xss_im&x/wb 
NE24&X  a  B22&x*Xm_im&x/wb 
NE42&X  a  -B41  &x*X88_im&x/Wb 
NE45&X  « -B41  &x*Xm_im&x/wt) 
NE51&X-  BS2&x*XssJm&x/wb 
NE54&X-  B52&x*Xm_i  m&x/wb 
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ND12&X  -  -8144x*XmJn>&x/#b 
N015&X  -  -8144x*XrrJmA»*to 
N021&X-  B254x*Xm_im&x*b 
N024&X-  B254x*XiT_im&x/Wb 
ND426x  - -B44&x*Xm Jm&x/wb 
ND46to--B444x*Xir_lm&*^ 
N051&X-  B56Ax*Xm_inr>&x^Nb 
N054&X-  B55&x*XrrJm&rfwb 

MACRO  END  krfimcoM 

imooaf  ("1”) 
imooaf  C2”) 
imcoef  (*3*) 

HimtII — Initial  conditions - 


1— murcM  and  loads 

iqsic 

-0.0 

-0.0 

iosic 

-0.0 

Hcqic 

-0.0 

jftflC 

*  1/Xmd 

Medic 

-0.0 

wric 

-  376.991 

kUpic 

-iqsic 

idpic 

-kteic 

topic 

*  KMC 

vqsk 

-1.0 

vdsic 

-0.0 

voeic 

-0.0 

mine 

-0.0 

! — exciter - 

vfdic 

-1. 

vreic 

»  vWc  ♦  1*exp(.3*Vfdfc) 

volic 

>0.0113 

vfbic 

-0.00008 

! — prime  mover — 

Tile 

-0.0 

wftic 

-  (THc/Clgt)  +  C2gt 

wfvic 

-  wftic 

werr3ic 

-wfvic- WT1  Os 

END  I  of  INITIAL 
DYNAMIC 

DERIVATIVE 

CONSTANT  vref  -1. 

CONSTANT  wref  -  376.991 
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!— Load  Parameters - 

CONSTANT  mot_on  1-0.0 ,  mot_on2-0.0 ,  mot_on3-0.0 
CONSTANT  cbl-0.0 

! — Square  lew  pump  loads - 

Ml  -  (wr1)~2 

02  -  (w»2)**2 

03  -  (wr3)**2 

!— R-L  load  parameters - 

CONSTANT  rl  -  5„  XI  -  2. 

! — Derive  quantities  for  output - 

zero  *0.0*1 

petec  ■  3*vas  *  ias 

pmech  *  O.SNwffVwb 

pdevt  *  0.5*wr*Te 

delta  ■  ias  *  IqId  *  id 

deNd  *  ids  -  kip  -  Id 

iamag  *  sqrt(iqs**2  +  kto**2  +  loe**2) 

IF  (ABS(iqs)  LT.  0.0001)  THEN 
iaphs-0.0 

ELSE 

taphs  «  atan(ida/iqs)*1  SO.Q^pi 

END  IF 

vrip  *  ABS(.07*coe(wr*t))  -  .035  !to  simulate  rectifier 
vamag  *  aqrt((vqs**2  +  vda**2  +  voe**2))  +  vrip 
frq  *  (wr  -  wbywb 

I— Convert  to  ind  motor  base - 

tl  *  te Jml  *smpb/impb1 

t2  *  te Jm2*8mpfaflmpb2 

13  *  te_im3*smpb/lmpb3 

wrl  *  wr_im1/wb 

wt2  *  wr_im2/wb 

wr3  » wrjm3/wb 

!!!!!!!!!— exlter  model - 

verr  *  vref  -  vamag  *  vfb 
vred  *  (verr*Kae  -  vreyTae 

vre  -  UMINT(vred,vreic,0.0,8.5) 

sat  *  vre  -  .1  *exp(.3*vfd) 

vfdd  *  (sat  •  vfd*Kee)/Tee 

vW  -  INTEG(vfdd,vfd»c) 

void  *  (vre*Kfe  -  voiyTfel 

vol  -  INTEG(vo1d,  voile) 
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vfbd  -  (void  -  vt>yTf«2 
vfb  « INTEG(vtbd,  vfoic) 

!!!!!!!!!— prime  mover/govemor  model- . . — 

CONSTANT  Gp-1. 
warn  •  Gp*(wref  -  wrywb 
werr2  «  KcNverrtflc  •  KcNwdtob 
werr3  » INTEG(werr2,werr3ic) 
wmnt4  -w©fr3  +  Wt10« 

WTvd  >  (wecr4  -  wfvyTTv 
Wfv  -  INTEG(Wfvd,WMc) 

wnd  -(wtv-wnym 
Wft  -  INTEGKWIW.WIWc) 

WR2  « (Wit  -  C2gt)*C1  gt 
Ti  -  WI12  +  CgngT warn 

I— derive  angles - 

third  «wr 

thtr  -  INTEG(thtrd,thtric) 

del  -atan(vdefcqsri80.Qfoi 

!!!!!!!!! — state  equations  for  synchronous  machine — 

iqsd  .  LI  1*lqs  +  N12*wr*lde  +  L14*kq  +  N15Me**d  & 

+  N16*wt*8cd  +  Virvqs 
iqs  •  INTEQ(lqsd,lqsic) 

Usd  ■  N21*Vff*iqs  +  L22*tds  ♦  N24*wr*kq  +  L25*Rd  & 

+  L26*flcd  +  V22*vds  +  V2SMd 
ids  -  INTEGfldad.idsic) 

toed  «  L33*ios  V33*Vos 

tos  « INTEGKtoed.toeto) 

ikqd  ■  L41*iqs  +  N42*wr*kJs  +  L44*Hcq  +  N45*wr*Rd  & 

+  N46*wr*Rcd  + V41*vqs 
ikq  .  INTEGflkqd.leqic) 

ifdd  a  N51*wr*iqs  +  L52*ids  +  N54*wr*Ncq  *  L55*td  & 

4-  LS6*Hcd  4-  V52*vds  4-  V55*Vfd 
Rd  -  INTEQ(lfdd.lMfo) 

Redd  -  N61*wr*iqs  4-  L62*ids  4-  N64Nw*toq  4-  L65*lfd  & 

4-  L68*Hed  4-  V62*vds  4-  V65*vfd 
Red  a  INTEGflkdd.kdlc) 

1— GENERATOR  oloctricaJ  torque  equation  In  teems  of  currents 
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Te  .  (Xnd*(^d»f«+*cd)*ltpHXmq*  (-iqe+fcqjlde) 

!— final  state  equation 

wrd  -wb*(Ti-Te)/2.(VH 

wr  ■  INTEG(wrd,wric) 

! — Speed  of  ynwtof  rotor  determines  iMricil  fnqutncy 
we  *  wr 

!!!!!!!!!— state  equation  macro  for  the  induction  motor  loads 
MACRO  imetateqn  (x) 

iqd_&x  -  LM11&x*iqJm&x  +  NE12&x*we*ldJm&x  +  & 

ND1 2&x*wdjm&x*id_lm&x  +  LM14&x*iqrJm&x  & 

+  NE1 5&x*we*1drjm&x  +  ND15&x*wdJm&x*ldrJm&x  & 
+  B1 1  &x*(vqe+vql)*mot_on&x 
iq_im&x  •  INTEQ(iqd_&x,0.0) 

idd.&x  ■  NE21  &x*we*iqjm&x  4-  ND21&x%rtJm&x*iqJm&x  & 

+  LM22&x*idJm&x  +  NE24&X*we*iqr Jm&x  +  & 
ND24&x*wdJm&x*iqrJm&x  4-  LM25&x*idrJm&x  & 
B22&x*(vd»4-vtM)*mot_on&x 
id_im&x  » INTEG(idd_&x,0.0) 

iod_&x  «  LM33&x*toJm&x  4-  B33&X*(voe4-vol)*moLon&x 
k>_im&x  -  INTEG(iod_&x,0.0) 

iqnj_&x  -  LM41  &x*iqjm&x  4-  NE42&x*we*kUm&x  +  & 

ND42&x*wd_im&x*idJm&x  4-  LM44&x*iqr Jm&x  4-  & 
NE45&x*we*idrJm&x  4-  ND45&x%KUm&x*idr_im&x  & 

4-  B41  &x*(vqs4-vqB)*mot_on4x 
iqrjm&x  -  INTEQ(iqrd_&xtO.O) 

idrd_&x  «  NE51  &x*we*iqjm&x  *  ND51 4x*wd_im&x*iq_im&x  & 

4-  LM52&x*kj_im&x  4-  N  E54&x*we*iqr Jm&x  4-  & 
ND54&x*wdJm&x*iqrJrn&x  4-  LM55&x*kJr Jm&x  & 

4-  B52Ax*(vds4-vci)*mot_on&x 
idr_im&x  « INTEG<idrd_&x,0.0) 


iord_&x  .  LM66&x*ior_im&x 
tor  Jm&x  -  INTEQ(tord_&x,0.0) 

Tejm&x  ■  Xmjm&x*(iq_im&x*kirjm&x  -  W_&n&x*iqr_im&x) 

wrd  Jm&x  «  wb*(T e_im&x  -  (Ti&x^mpb&x/8mpb)V2.0/H_im&x 
wrjm&x  -  INTEQ(widJm&x,0.0) 
wd_im&x  »  we  •  wr_im&x 


97 


MACRO  END  lofimotatoqn 
tilMHi— cal  induction  motor  modal  macro 


imatalaqnO*) 

imatatagn  (*2*) 
imatatagn  (*3*) 

!!!!!!!!! — sum  of  motor  cunents  and  current  derrivatNes  for 
t— -conatraint  aquation 

iqtd  -  (iqd_1  ♦  iqdL2  ♦  kjd_3) 
kjl  » (iqjml  +  iq_im2  +  iq_im3) 

idU  -  (Wd_1  +  Wd_2  +  WCL3) 

U  - (kLiml  ♦  W_kn2  +  WJm3) 

iokf  « (iod_1  +  tod_2  +  tod_3) 
iot  « (tojml  +  tojm2  +  to_im3 ) 

!!!!!!!!! — state  equations  for  the  paraM  R-L  load - 

kftxj  « (-wb*rVXI)*iqlp  -  wTk%>  +  (\wtyXI)*(vq8+vql)*cb1 
k*P  =  INTEG(iqlpd,kflpic) 

idlpd  .  wr*iqip  •  (wb*rVXI)*k«p  +  (wtVXI)*(vd8+vdl)*cb1 
idp  =  INTEG(«pd,idpfc) 

tofcd  « (-wb*rVXJ)*totp  +  (wtVXI)*voe*cb1 
to *>  -  INTEG(totpd,toipfc) 

!!!!!!!!! — DASSL  bus  volage  model  based  on  impact  relation - 

f»»q  -iqad  +  iqs-iqlpd-kfr-iqkl-iql 
vqs  « IMPLC(re8iq,vq8fc) 

resid  «  tied  +  kte  -  kJpd  -  •  kfld  •  Ml 

vds  >  IMPLC(reeto,vdsto) 

resio  «  toed  +  toe  -  toipd  •  tolp  -  told  -  toi 
vos  m  IMPLC(resto,vosto) 

IllUilll — Hne  toss  model - 

CONSTANT  rll » .005 
CONSTANT  xll  -  .001 

vq*  «<-  rTiqe  -  xffqed/tob  •  wrxTiqsAvb) 
vdfl  *(*  rtrkte  •  xtndacVwb  +  wfxTiqs/wb) 
vd  -(-  rtTtos  •  xtrtostVwb) 
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!!!!!!!!!— convert  currents  to  a-b-c  reference  for  output - 

iae  -  iqs*coe(thtr)  +  kte*sin(thtr)  +tos 

ibs  -  iqs*coe(thtr-2.0*pt/3.0)  +  Ws*stn(thtr-2.0*pl/3.0)  +ios 

ics  *  iqs*cos(thtr+2.0*pi/3.0)  +  ids*$in(thtr+2.0*ptt3.0)  +  ios 

ial  =  jqTcoe(thtr)  +  kJI*sin(thtr)  +  iol 

Ibl  >  iqi*cos(thtr-2.0*pi/3.0)  4-  k**sin(thtr-2.0*pi/3.0)  +  iol 

id  « iql*co8(thtr+2.0*pi/3.0)  +  idl*sin(tMr+2.0*pV3.0)  +  iol 

!!!!!!!!! — convert  voltages  to  a-b-c  reference  for  output - 

vas  =  vqs*co8(thtr)  +  vds*sin(thtr)  +vos 

vbs  *  vqs*cos(thtr-2.0*pi/3.0)  +  vds*sin(thtr-2.0*pi/3.0)  +  vos 

vcs  =  vqs*cos(thtr+2.0*pi/3.0)  +  vds*sin(thtr+2.0*pi/3.0)  +  vos 

END  !  of  DERIVATIVE 
CONSTANT  tstop  =  6.25 
TERMT(t  GE.  tstop) 

END  I  of  DYNAMIC 
END  !  of  PROGRAM 


!!!!!  --COMMAND  RLE  FOR  DASSL  SYSTEM  MODEL—  lllllll 

prepare  t,jqs,iqjm1  ,vqs,vds,vfd,Te,wr,del,ias,vfb,vre,verr,werTl  ,ti,frq  & 
deliq, delid, iqsd.iqid, zero, vamag.iamag.iaphs, vas, tl  & 
iqjm2,iqjm3,iqlp, t2.wri  ,wr2,pelec,pmech,vql,vdl 
set  calplt=.f.,strplt=.t.,alcplt=.f. 
set  nrwitg=.t. 

proced  mot 
start 

s  mot_on1  =1 .  mot_on2=1 .  mot_on3=1 .  tstop=9.0 

cont 

show2 


proced  mot2 

s  mot_on1=1.  mot_on2=1.  mot_on3=1 .  tstop=3.0 
start 


proced  stpid 
start 

s  mot_on1  =1 .  mot_on2=1 .  mot_on3=1 .  tstop=9.0 
cont 
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•  rt-2.  xqW2.  xdW2.  xot-2.  cbl-1.  tstop-12. 
cont 
ahow2 
and 

pfoced  short 

s  rl-.OOOl  xdU  OOOl  xqt.0001  xoi-.OOOl  ct>1«1.  tstop-6.5 
coni 

s  rUI 5.  xdMO.  xqklO.  xoMO.  cbl-O.  tstop-8.0 
cont 

plot  las  deliq  vas  vamagfrk>-6.2/xhM8top 
pause 

ptotfrqti  vfd 
end 

procedshow 

plot  wr2  t2/to*-5/hi=5  wrl  t1/k>»-5/hi^/xk>o6.Cyxhl*tstop 
pause 

plot  iaa/los-1 ,5/hi»1 .5  deliq  vas  vamag 
pause 

plot  frq/to=-.01/hi-.01  ti/k)«-.2/hU.4  vfd/lo»0.0 
end 

proced  savplt 

s  xtteph. 16667  yttepl*.2  grdspt.t.  satspl-.t. 

sdevptt=4  ptt=11 

plot  wr2/hi=1 7xk>»6.0/xhUt8top 

s  plt*12 

sytispU.16667 

plot  t2/to=-6/hU6 

splt=13 

s  ytispl=.2 

plot  wrl /hM. 

s  plt=14 

plot  t1/k>=-5/hla5 
splt-15 

s  ytwpfcr. 16667  yinspUl  .6667 
plot  ias/k>=-1 ,25/hUI .25 
splt>16 

s  ytispU.25  yinspU2. 
plot  vas/k)»27hls2. 
splt-17 

plot  fiq/k>=-.01/hU.01 

splt»18 

sytispU.2 

plot  tiflo--.5/hU.5 

splt»19 

plot  vfd/kuO.O 

spft-20 

s  ytispU.2  yinspUl  .4 
plot  vamajytiUI  .4/1o*0.0 
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end 


syinspi-2. 


proced  show2 

8  xttepU. 16667  ytiapU.2  grdapU.t.  satapU.t. 
8devplt-6 

plot  wr2/hU1  Aj-OTxio-e.O/xhWstop 
pause 

syttepW. 16667 
plot  t2/to~-6/hi-6 
pause 
s  yttepls.2 
ptot  wr1/hi=1. 
pause 

plot  t1/to*-5/hU5 
pause 

s  yttepU.25  yinspW2.5 
plot  ia8/k>*-1 .25/Wal  .25 
pause 

s  ytispl«.25  yinspl»2. 
plot  vas/to=-27hi=2. 
pause 

plot  frcyk)*-.01/hi=.01 

pause 

s  ytispl=.2 

plot  ti/lo=-.5/hl=.5 

pause 

plot  vfd/k>=0.0 
s  ytispl=.2  yinspM  .4 
pause 

plot  vamag/hi=1  4/k>=0.0 
s  yinspt=2. 
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